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INTRODUCTION 


The determination of the nuclear radiation environment at some 
point in a complex system is a difficult and often costly procedure. 
However, because of the harmful effects that nuclear radiation has on 
both organic systems such as man and inorganic systems such as elec- 
tronic components and structural material, it is necessary to deter- 
mine the environment in all systems that are exposed to any significant 
amount of nuclear radiation. This is particularly true for certain sys- 
tems which have been of great concern to the National Aeronautics and 
Space Administration (NASA) over the last 15 years. Examples in- 
clude nuclear propulsion systems such as NEF.VA (Nuclear Engine 
Rocket Vehicle Applications) and RIFT (Reactor In Flight Test). The 
SNAP (Systems for Nuclear Auxiliary Power) program resulted in the 
design and construction of both nuclear and radioisotopic power plants. 

51 AP type power systems have been used for powering experimental 
pacaugcc on the lunar surface (SNAP -27) and for the Nimbus satellite 
(SNAP-19). 

At the present time, NASA's interest in nuclear systems is 
much less than it has been in the past. However, it is likely that in 
the future many missions, especially those involving manned inter- 
planetary exploration, will require nuclear systems. For example, 
a manned Mars mission would be very difficult to accomplish with pre- 
sent chemical propulsion capabilities, but almost all of the technology 
necessary for a Mars mission using a N 1RVA or NERVA type propulsion 
system is available. Using nuclear systems for iarth escape, braking, 
and return make the Mars mission a possibility n the 1980’ s. One 
aspect of the technology that requires further refinement and which 
also impacts all phases of the mission is the radiation environment 
produced by the nuclear propulsion system. This report discusses a 
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study made of a technique to improve the calculation of the radiation 
environment. The technique is called the iterative forward -adjoint 
Monte Carlo method. 

For a large class of radiation transport problems, only the 
Monte Carlo method has proven to be a useful tool for effecting solutions. 
This is primarily due to the fact that three -dimensional, time dependent 
problems can be solved by the Monte Carlo method, while other methods 
are usually restricted to only two dimensions. Monte Carlo suffers 
from the disadvantage that as a probabilistic method, considerable 
computer time is required for statistically meaningful results. This 
limitation can be alleviated by the use of variance reduction tecnniques. 
The research that has been performed in this study has identified 
variance reduction techniques based on both the forward or normal 
Monte Carlo solution and the adjoint Monte Carlo solution which can 
be used to reduce the variance for a given amount of computation time. 

In particular, the techniques utilize the fact that the adjoint fluence 
determined in an adjoint Monte Carlo calculation can be used to calculate 
altered sampling distributions (for the source parameters, path length, 
and past -collision parameters) to be used in the forward calculation. 

The source term used in the adjoint calculation must be the detector 
response function in the forward case. Random sampling from these 
altered distributions, instead of the natural distribution, results in a 
reduced variance of the effects of interest because the most "important" 
part of the distribution is sampled most often. The forward fluence 
determined in a forward calculation can likewise be used to alter the 
distribution function for the adjoint calculation. The validity of these 
techniques has been discussed by many researchers (e. g. , see 
References 1-8). The earliest reference is contained in Herman Kahn's 
Application of Monte Carlo (1), first published in 1954. He also 
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suggested the use of the iterative forward-adjoint Monte Carlo (called 
Method IV), stating "As far as the author knows, Method IV has never 
been used in a systematic fashion." By iteratively performing the 
forward and adjoint calculations, updating the altered sampling distri- 
butions between each calculation, the /ariance should reduce more 
rapidly with each iteration, since be cer values of the fluence (forward 
and adjoint) are being determined. 

The principal effort of this study has been to develop these 
variance reduction techniques which can be applied to computation 
methods, and a computer test bed written which will perform the 
appropriate calculations for testing these different techniques of using 
the forward or adjoint fluence to calculate the altered sampling 
distributions. 

The general outline of an iterative forward -adjoint Monte Carlo 
calculation is given below and diagrammed in Figure 1. It is assumed 
that the particular problem of interest has been identified and the 
necessary data obtained and p~epared in a form acceptable to the 
computer program which employes the iterative forward -adjoint 
variance reduction techniques. A discussion of such a program, 
called IFAM for Iterative Forward-Adjoint MORSE, is given in 
Chapter 3. The initial step in the calculation is the processing of input 
data. This step is performed for both forward and adjoint data before 
any of the random walk calculations begin. Data which is not in use 
(e.g. , adjoint data during the forward mode random walk) is stored 
on bulk storage or data files units such as disk files. The input data 
must also specify the initial mode - forward or adjoint - to be executed, 
since this is left to the judgment of the user. 

After the processing of the input data, the data for the initial 
mode is retrieved from the proper data file, data initialization performed 
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^igure 1 . Iterative Forward -Adjoint Monte Carlo 
Flow Diagram 
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and the data put into core. The initial mode, either forward or adjoint, 
random walk is begun for the specified number of batches and histories 
per batch. Ekiring this calculation, the program executes only the 
initial mode type of random walks, but it employs input biasing param- 
eters, if desired, to alter the source particle energy and direction, 
the path length, and the energy downscattering distribution functions. 

At the same time, data is stored for estimating the energy and angular 
dependent flu> .ce in each importance region. Estimates are also 
made of some effect of interest, 3uch as the dose at the detector. This 
process is continued until all histories have been completed for the 
initial mode, at which time the output data is written, including the 
estimate of the effect of interest and the energy and region dependent 
fluer.ee estimate. The data used or generated during this calculation 
is stored back on the data file uiiits for processing and use in the next 
iteration. 

The calculations performed in the second half of this first 
iteration are very similar to those performed for the initial mode ex- 
cept that the other mode type random walk will be executed and the 
distribution functions will be altered with both the input biasing functions 
and the energy, angular and region dependent fluence from the initial 
mode. The validity of using the initial mode fluence for altering the 
distribution functions is discussed in detail in Chapter 2, and the form 
employed in this study is explained in Chapter 4. Retrieval of the 
opposite mode data from the data file; initialization and generation of 
the importance function is required before the random walk calculations 
can occur. At the completion *? these calculations, the same results 
as in the initial mode are output both to the printer and to the data file 
for this mode. This completes the first iteration which consists of 
random walk calculations in both the forward and adjoint modes. Note 
that at the end of each mode, an estimate of the effect of interest was 
obtained. 
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The second and succeeding iterations are identical in format 
to the first iteration: data is retrieved for the initial mode data file, 
initialization and importance function calculations are performed, thtn 
the effect of interest and the fluence is estimated during the random 
walk calculations followed by output and data storage for the initial 
mode The same steps are then performed for the opposite mode. 
However, with each iteration, the fluence estimates approach the 
actual value, thus improving the altered distribution functions and 
hence reducing the variance of the effect of interest estimator with 
each iteration. An evaluation of the iterative forward-adjoint technique, 
including some general guidelines for employing the technique, is given 
in Chapter 5. 
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2 . 


THEORETICAL ASPECTS 


The fundamental equat^n describing the -ransport of radia- 
tion through matter is the Boltzmann transport equation. This linear 
integro-differential equation is essentially a particle balance equation 
in phase space. In Section 2. 1, each term of this equation is defined 
by discussing both the physical significant and the mathematical form 
of each term. The physical assumptions upon which the transport 
equation is based and which restrict its application are also presented. 
Whereas the Boltzmann transport equation provides the basis lor 
numerical solution methods such as discrete ordinates and the moments 
method, an integral transport equation (i.e. Fredholm equation of the 
second hind) is a better form for stochastic or Monte Carlo method 
solutions. However, the integro-differential form seems to give a better 
physical insight into the radiation transport problem, and for that rea- 
son the Boltzmann equation was chosen for the introductory discussion 
of the iterative forward -adjoint Monte Carlo method. 

Section 2. 2 presents a derivation of the integral transport 
equation from basic definitions. The concept of the emergent particle 
density (or the densitv of particles leaving sources and collisions) is 
also introducted. This concept is important to the Monte Carlo method 
because the emergent particle density equation is the one simulated in 
most Monte Carlo computer programs. The format of the presentation 
in Section 2. 2 and Section 2. 3 follows that of the report by Goertzel and 
Kalos (Ref. 2), which is an amplification of the fundamental work done 
by Kahn (Ref. 1). Section 2. 3 illustrates the Monte Carlo method by 
an example transport game, and then uses this game to produce the 
’'perfect game" by the introduction of the adjoint or importance function. 
The implications of these results are then discussed relative to the 
iterative forward-adjoint Monte Carlo method. 
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Since the computer test bed used to validate the iterative for- 
ward-adjoint variance reduction techniques is a multigroup Monte 
Carlo code, Section 2.4 contains a discussion of the multigroup in- 
tegral transport equations in both the forward and adjoint forms. 

These equations are used to define the proper relationships between 
the forward and adjoint forms. 

2. 1 Boltzmann Transport Equation 

The Boltzmann transport equation (sometimes called the linear 
transport equation, Boltzmann equation, or transport equation) describes 
the distribution of particles such as neutrons or gamma- rays in phase 
space, P. Phase space will be represented by: 

P = (x, E,t) = (it, E, co,t) 

where x = spatial variable (e. g. x = (x,y, z) ) (vector) 

E = directed energy variable (E = (E, w) ) (vector) 

E = energy variable (scalar) 

co = direction or angular variable (to = (0,<p))(vector) 
t = time variable (scalar) 

The Boltzmann equation is usually written in the following form: 

[l It * 0 + v *6,E,t) + 2^(E)*(x,E,t) 

J* J 2(x,E') * (x, E t) f(E,i|E',u')dE'du»' 
s 

+ S(x, E, t)] dicdE (2-2) 

= differential volume element (e. g. dxdydz). (scalar) 

= differential energy variable (scalar) 

= differential angular variable (e. g. sin 0ded«p) (scalar) 


where: 

dx 

dE 

dw 


( 2 - 1 ) 
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dE 


= dEdto 


$(x, E, OdxdEdco 


^(x, 12) 

E s (x, E) 

f (E,w|E',w')dEcJw 


S(x, E, t)dxdE 


= particle flux at time t in the differential 

volume dx about x with energies in dE about 

E and with directions in dw about to. The 

2 

units of <J> are particles (or photons)/cm /sec/ 
steradians/eV. , and the energy E corresponds 
to the speed v. 

= macroscopic total cross section in cm * at 
position x and ener,:v E. 

= macroscopic scattering cross section in cm 
at position x and energy E. 

= the conditional probability of a particle 
scattering from an initial energy E' and 
direc.ion uj' into the energy interval dE about E 
and into the £.jhd angle dco about the direction to. 
= the radiation at time t in the differential volume 
dx a!, jilt x with energies in dE about E and with 
directions in dto and to appearing from sources 
other than scattering events (e. g. fission 
neutrons, radioactive decay gamma-rays). 


The first term on the left side of the Boltzmann equation 
represents the time rate of change of the radiation in the differential 
phase cell, dxdE, a time t. This time rate of change can be caused by 
the following four occurrences: 

1) Leakage [to * v^x^E, t)dxdE] : Net leakage of particles 
with directions in dto about to and with energies in dE 
about E from the differential volume, dx about x, at 
time t, 


(scalar) 

(sc alar) 
(scalar) 
(scalar) 

(scalar) 

(scalar) 
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2) Interactions [2^(x, E)<t(x, E, t)dxdE] : Any interaction of the 
particles with the medium at time t in dx about x which 
removes the particle from cither dE about F, or do; about oo. 


3) Inscattering / j* 


r 


£ s (x, E)«^x, E, t) f(E, u.|E', w')dE 'dw ' ) dEdw*| 


A gain in radiation in dx about x at time t due to any scatter- 
ing events which redirect the particles into do; about a> and 
the energy into dE about E, 


4) Sources [S(x, E, t)dxdE] : Particles born in dx about x with 
the proper energies and directions as defined above. 


These four terms constitute the losses (leakage and interactions) 
and gains (inscattering and sources) which determine the time rate of 
change. Thus, the Boltzmann equation is essentially a particle balan-’p 
equation, that is: 

Change = Gains - Losses. 

It should be noted that equation 2-2 is not the only form of the Boltzmann 
equation. Often the balance is performed on the phase space particle 
density instead of the particle flux. The independent energy variable E 
can be replaced by the speed, v. Also, the inscattering term, which 
includes any scattering event in which the number of particles is not 
changed, can be replaced by a more general term. This term can include 
production events such as fission in neutron transport and pair production 
in gamma-ray transport. When the inscattering term is thus modified, 
the source term must be modified corresponding. 

The justification of the linear Boltzmann transport equation 
as given in Equation 2-2 requires that the following assumptions be 
made (also see Reference 9): 

1) Statistical fluctuation are neglected 
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2 ) Collision times are negligible 

3) Particle correlations are neglected 

4) The wave nature and spin of the particle is mglected 
(except in the computation of the cross sections and 
scattering kernel, f(E, o> |E u/) ) 

5 ) Interactions between particles are negligible 

6) The medium in which the particles is being transported 
is not affected by the presence of the particles (i. e. the 
cross sections and scattering kernel are assumed not to 
depend functionally on 4>(x, E, t) ) 

7) Cross sections are assumed independent of co 

8) The scattering kernel, f(E, co| E ', depends only on the 
angle between and u ' (or on w * a; ') and on E*. 

Assumptions I) through 4) are required since the Boltzmann equation is 
derived in the continuous phase space P based on classical physics. This 
means that the independent variables, x, E, oo, and t, are assumed to be 
continuous, not discrete (however, most transport methods use discrete 
approximations for some or all of these variables). Assumptions 5) and 
6) were necessary to keep the Boltzmann equation linear. Thus Equation 
2-2 can be used for neutron, gamma-ray, and charged particle transport, 
but would fail to adequately describe low energy X-ray transport because 
of stimulated emission and the extreme temperature sensitivity of cress 
sections. Finally, assumptions 7) and 8) avoid the necessity of having 
a preferred coordinate system. These assumptions would not be valid 
for the Bragg (coherent) scattering of low energy neutrons by crystals. 

The integral form of the transport equation can be derived 
from the Boltzmann equation by defining a new quantity of interest, the 
emergent particle density (Ref. 10). A derivation of the integral transport 
equation in multigroup form has also been derived for the MORSE code 
(Ref. 11) from the Boltzmann transport equation (Eq. 2-2). These 
derivations will be discussed in Section 2. 4. However, a derivation 


2-5 



of the integral transport equation from more fundamental considera- 
tions is given below. 

2. 2 Integral Transport Equation 

The following discussion of *he integral transport equation 
is valid for both neutrons or gamma -rays, provided that the assump- 
tions imposed oi> the linear Boltzmann equation are valid. Each par- 
ticle will be described by its position in phase space (P), consisting 
of the spatial coordinates, x, energy, E, and direction of motion, Z > . 

Only time independent cases will be considered, but addition of time 
dependence is straightforward. Hence, P can be represented as 
shown below: 

P = (x, E) = (x, E, w) . (2-3) 

The motions of the particles will be described by three terms: 

• Flux - $ (P) 

• Collision density - <A(F) 

• Density of particles leaving collisions (emergent particle 
density) - x (P). 

Note that all three of these terms are functions of phase space, P, and 
not just x or (x, E) . As usual, the flux and collision density are related 
by the total cross section, 2^(x, E), which is assumed to be independent 
of to, so that 

<Mx, E,u>) = $(x, E) t(x,E,w) . (2-4) 

In order to derive the transport equatio n, it is convenient to 
consider the number of collisions that a given particle has undergone. 
Therefore, the following definitions are in order: 

• *n ^ “ ^ ux a * ^ °* P ar ticles that have undergone n - 1 
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collisions (hence are entering their nth collision) 

• i// n (P) - collision density at P of particles that have under- 

gone n - 1 collisions fas for (P)] 

• X n (P) - density of particles at P which have jus; had 

their nth collision 

• x Q (P) - density of particles at P which are emitted from 

the source (consider the source as the zeroth 
collision) 

• C (E|E 7 ; x) dE - probability that a particle having a colli- 

sion at the spatial point x and directed 
energy E 7 will emerge in dE about E 
[C (E| E 7 ; x) is called the collision kernel] 

• T (x|x 7 ; E) dx - probability that a particle leaving a colli- 

sion at x 7 with directed energy E will have 
its next collision in d>; about x [T ix|x'; E) 
is called the transport kernel ]. 

From the above definition, the following relationship can be deduced; 


* (P) = X ( p ) 

n=l 1 


ao 

v 


<MP) = X W 

n=l n 

X(P) *<W 


00 

x 0 (p) + y x_ (p) 


k * 


X 0 (P) = S (P) (the source at P) 


(2-5) 

( 2 -«) 

(2-7) 

( 2 - 8 ' 


X n (x, E) = / C(E|5';x) \ (x,E') dE" n = 1,2,... (2-9) 

«J/ Ax,E) = [ T(x|x'; E) X (x(E)dx' n = 0,1,2,... (2-10) 

n+l J n 
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Note that in Equations 2-5 and 2-6, the summation begins at n ~ 1, 
but in Equation 2-7, it begins at n = 0. The reason for this is that 
$ (P) and (P) represent the particle immediately before entering a 
collision, and since n = o denotes the source, $ o (P) or ^ (P) are 
physically meaningless. 

The physical significance of the transport kernel can be seen 
from the fact that T (x|x'; E) is just X^(x, E) exp [I^(x, E) • jx - x' |] 
if the material property from x' to x does not change and if the vector 
frem x' to x is parallel to tu, where E = (E, ta. If the vectors are not 
parallel, then T (xlx' ; E) is identically zero. Thus the significance 
of Equation 2-10 is that the collision density at P = (x, E) for particles 

entering their (n+l)th collision is just the integral over all spatial points 
of the density of particles leaving x' with directed energy E times the 
probability that the particle will reach x without having a collision. The 
physical significance of Equation 2-9 is very similar except that the particle 
entering into the nth collision dots not change spatial coordinates, but 
changes its direction from co’to £ and its energy from E’ to E. Thus 
Y n (x, E ) is just the integral of the collision density times the collision 
kernel over all directions and all energies. The exact form of the collision 
kernel is much more difficult to describe than the transport kernel since 
it must include all possible types of collision events, such as elastic and 
inelastic scattering, absorptions, and fissions, whenever applicable. 

By substituting Equation 2-9 into Equation 2-10, where x is 
replaced by x’ , then y n (x’ , E) can be eliminated: 

^ nfl (x,E) = J J T (x|x'; E) C(E|E';i') *„(*',£') dE' d k' 

n - 1,2,... (2-11) 


For n = o, we use the fact that ^(P) = S(P), so that 
0 1 (x,E)= J T(xix' ; E) S(x, E) dx’ . 
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Note that the order in which the collision and transport kernels appear 
is very important, s. nce the collision density must first be transformed 
by the collision kernel before the transport kernel can operate on it. 
Substituting Equation 2-10 into Equation 2-9 for the (n+l)th collision yields: 

*n + 1^* ’ E ^ = 7 J c (^|E ;x) T(x|x'; E ) X" n (£', E'ldx' dE' 


n - 0, 1, 2 

In Equation 2-13, the order of the collision and transport kernels are 
reversed to that Equation 2-11, as one would expect from physical 
considerations. 


(2-13) 


The relationships between ^ (P) and Y (P) can be derived by 
summing Equation 2 3 over n = 1,2,...! 

CD oo 

I x n (i, E) = g^ j C(ElE';i)^ n (x, E'l dE' 


(2-14) 


or: 


°o oo 

(i,E) -y (i,E) / £c(EIE';xW ( -j) ^ (2-15' 

n=o n n=l 


By Equations 2-6, 2-7, and 2-8, Equation 2-15 can be written 
X(P)=j C(Ej2*;x)tJ> (x, E ) dE' + S(P) 


(2-16) 


Likewise, by summing Equation 2-10 over n = 0, 1, 2, ... ! 




or: 


*Hx, E) = J* T(x)x y ;E) X (x* £)dx' . 


(2-17) 


(2-18) 
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Substituting Equation 2-18 into Equation 2-16 will clear Equation 2-16 of 
4»(x, E') as follows: 


X (x, E) = j" { C(E|E / ; x) T (x|x'; E') V (i; E') d x’ dE’ 


+ S(x, E) 


(2-19) 


Likewise, X(x'. E) can be cleared from Equation 2-18 by substituting 
Equation 2-16: 


4> (x, E) T(xlx';E) C(EjE';x )«Mx',E') dE'dx' 

+ I T(x|x‘; E)S(x', E) dx' . 


(2-20) 


Equation 2-20 can be written in the form: 


‘J'(P) K(P(P')0(P') dP* + Q(P) 


( 2 - 21 ) 


where: 


K(p|p ) = T(x|x':E)C(E!E';x') (2-22) 

Q(P) = J T(i|x';E) S(x',E) dx' (2-23) 

dP' = dE' dx' . ,2-24) 

Equation 2-19 and 2-21 are two forms of the integral trans- 
ports equation. The equation most often simulated in forward Monte 
Carlo calculations is 2-19. However, due to the fact that most people 
have abetter understanding of the collision density equation (2-21), 
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this equation will be used in the next section to introduce the concept of the 
adjoint to an integral equation, to define the term, adjoint, both mathematically 
and in a physical sense, and to illustrate how the adjoint equation solution 
can be used to improve the forward Monte Carlo estimate of an effect of 
interest. 

2. 3 Monte Carlo Transport Gam es 

Suppose one want to evaluate an effect of interest defined by: 

<F> =yV (P) f(P) dP, (2-25) 

where 0(P) is the collision density defined in Equation 2-17 and f(P) is a 
collision density response or payoff function (see Section 2.4.2 for a 
discussion of the form of the response functions). Since the calculation of 
the collision density by analytical techniques is impossible except for very 
simple problems, one must employ a numerical technique to evaluate or 
estimate the effect of interest. The Monte Carlo simulation of Equation 2-17, 
which consists of estimating the expected value of (2-25) by a series of 
random walks based on Equation 2-25 is one possible technique. This 
procedure, described in detail in Section 3. 2, consists of picking 
"particles” from the first collision source term, Q(P), (i.e. , selecting 
the spatial, energy, and directional parameters of the particle) and 
then determining the subsequent history of the particle from the kernel, 
K(P|P’). Certain restrictions are imposed due to the statistical nature 
of this procedure, for example, the first collision source must be a 
properly normalized probability function, either discrete or continuous. 

Also, Q(P) and K(P I P' ) must nonnegative. Thus, a source normalization 
factor is required for Q(P). Equation 2-23 indicates that the normali- 
zation factor is not the magnitude of the natural source, S(P), for a 
finite system. The distribution of the particle collisions (and any daughter 
particles) so transported is proportional to the magnitude of the collision 
density. Thus, the values of the response function, f(P), at each collision 


2-11 



is an unbiased estimator of Equation 2-25. An estimate, F, of <F> 
can be calculated for a finite number of particle histories, as shown below. 


_ i N C 3 , 

'< pi „> > (2 - 26 > 

j=l n=l 

where N is the number of histories generate! and C. is the number of 
collision which the j-th particle and its daughters experience. It lias been 
assumed that each history terminates is a finite number of collisions (i. e. , 
the process is subcritical and C. is finite). 


It is also possible to generate an unbiased estimator of <F> by 
choosing from a different or altered source and kernel (say Q and K), if the 
contribution from the response function at each collision is weighted by the 
product of the ratios of the true to altered source and kernel distribution. This 
weighting factor after m collision, v. T m , which multiplies the response value 
can be written as: 


Q(P.) m K(P P .) 

w_ ^ 1 n n n-1 

Q(P.) n=2 K(P P .) 

1 ' n n-1 

m 

■ W o (P l> ■ ", W l (P n| P n-l ) - (2 - 27 ' 

n-4 

Equation. 2-27 defines tne weighting functions which are used in Figure 2 to 
illustrate an altered transport game, in order to define the game properly 
the following requirements on our altered source distribution and kernel are 
imposed: 


Q(P)> 0 ) (2-28) 

/§( P)dP = l \ 
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SELECT k FROM p^P-.,) 
WHERE: 


Po< P n-1> " 1 Pl< P n-l‘ 

Pl<*Vl> " y* (P i P n-1> dP 

tV(P iVl ) - OlFk > 2 



Figure 2. An Altered Transport Game 
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and 


K (P | P’ ) > 0 


; 

/ K(P|P’)dP H Pj(F)<l ) 


(2-29) 


Equations 2-28 assure that the source term is a properly normalized 
probability f notion. Equations 2-29 restrict the transport problem to one 
in non -multiplying media, thus assuring the convergence of the transport 
in a finite number of steps if there exists any medium which has a positive 
absorption probability. This probability, p Q (P’), of entering a trapped 
state at P' is just: 


P 0 (P f ) = 1 - P^P’), (2-30) 

which is the absorption probability, and it is assumed that each random walk 
enters a trapped state. 

As a further modification to the transport game , instead of 

evaluating the response function, f(P’), at each collision event, score only 

when P' is a trapped state for the particle. Such a change in the scoring 

technique requires a change in the response function without changing the 

expected value. This can be done easily by considering that the density of 

trapped states or absorption event density, which will be denoted by 4 / . 

ft 

The absorption event density is related to the collision density by: 


<A a (P) = p o (P)<A(P) (2-31) 

where p Q (P) is the absorption probability. Thus, the effect of interest in 
Equation 2-25 could also have been calculated with the integrand, 

■A a (P) f(P)/p (P), and so the response function defined by f(P)/p (P) is also 
an unbiased estimator of < F> when only those collision resulting in an 
absorption contribute to the estimate of F. Combining this new response 
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function with our altered source and kernel produces: 

? -wi. ^ I < p i>^o< p i>- (2 - 32) 

. 1=1 

where w is the weighting function value calculated at the absorption event, 

3 . 

which occurred at r> ^ for the j-th particle. The flow diagram of this altered 
game is depicted in Figure 2. All terms are defined in the convention used 
above, except the particle superscript, j, has been dropped for clarity 
purposes. Now consider methods of reducing the variance of the estimate 
of < F> by utilizing the adjoint equation. 

2.3.1 The Adjoint Equation 

Although the transport game discussed above was based on 
altered probability distribution functions, no mention was made of the way in 
which these distribution functions are chosen (note that K(P| P')/p G ^P') 

is the distribution function, not K( T5 |P')). Obviously, one wants to alter 
the distribution functions in such a manner that scores (F) yield a 
better value of F and io do this with smaller sample size (i.e. , smaller 
N). Therefore, one seeks to choose "particles" from Q and points in 
phase space that lead to higher expected contributions to F. This means 
that one needs to have some sort of importance function, say I(P), which 
is at least approximately proportional to the expected value of F from a 
particle at P. The with I(P) we can define appropriate values of £$(P), 

P o (P), and K(PlP'). so that the variance of our calculation is minimized, 
where the variance is defined as 


V 





(2-33) 


To motivate our selection of the above values, consider the 
following discussion of the equation which is adjoint to Equation 2-21. 
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Since K(p|p ) i;3 the product of two real kernels, TCxix*"; E) and 
C(E|E ; x*), then K(P|P*) is also a real kernel and its adjoint is 
just K*(P jP < ) = K(P'IP) 

Therefore, any equation of the form 

J (P) = Jk(P' | P) J(P') dP' + i(P), 

assuming that J(P) and f(P) are quadratically integrable, is an adjoint 
of Eq. 2-21. Defining the operators: 

K*(P)=J K(P| P^(P') dP' 

. and K* J(P) = J K(P’| P) J(P) dP' 


the following rc.ationship is known to be true (Ref. 12): 

j"j(P) • K vHP)dP =J</' (P)-K* J(P) dP • ( 

The proof is straightforward, requiring only a reversal of the variable 
of integration and rearranging of terms inside the integrals. Multi- 
plying Equation 2-21 by J(P) and integrating over P and then multi- 
plying Equation 2-34 by </» (P) and integrating over P, produces another 
identity: 

j (P) f(P) dP = J J(P) Q(P) dP • (2- 

But from Equation 2-25 it is obvious that: 

<F> = j J(P) Q (P) dP • (2 

Thus, another method for determining < F > is to perform a transport 
game by Monte Carlo methods on the adjoint integral Equation (2-34) 


(2-34) 

(2-35) 

(2-36) 

2-37) 

•38) 

•3G) 
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aa^ estimate the average value of <F> by: 


< F > s 


\_ 

N 


N J 

EE 

j=l n=l 



(2-40) 


This method is analogus to the on? described earlier in which Equations 
2-21 and 2-26 were used. 

One other mathematical characteristic of the operator, K, and 
its adjoint, K*, is that since: 

K>A(P) - /t(x| x’ E)Jc(E | E’; x’)<Mx\ E') dE’ die’, ^-41) 

Then K* J(P) = fc(h}E.k) f T(x| x;£') J(x ', E') dx’ dE’. (2-42) 

This is due to the fact that the adjoint of the product of two operators 
is just the product of the adjoint to the operators taken in reverse order. 
This relationship will be used in Section 2. 4. 

Now consider the following physical identifications to those 
terms which are in Equations 2- 25 through 2-39 and were not previously 

identified. In Equation 2-25, the function f(P) shall be a response or 
detector function which produces some contribution to the effect of in- 
terest due to a particle that has a collision at P in phase sp* ce. The 
effect of interest < F >, is a physical quantity such as the dose, heat- 
ing rite, or fission per kilogram in Uranium-225 The distribution, 

Q(P), is just the so-called "first collision source, or the distribution 
of point at which particles leaving the real sources, defined by S(x, E) 
experiences their initial collision (see Equation 2-23). By choosing to 
use the response function, f(P), in the adjoint equation (2-34), it was 
possible to construct a new method of solving for the effect of interest. 
Mathematically, any suitable (i. e. , quadratically integrable) function 
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could have been used instead of f(P) in Equation 2-34, but f(P) was 
chosen because any other choice would not produce a second method 
of estimating < F>. 

The physical interpretation of our adjoint eq> ^tion, as given 
by Goertzel and Kalos (Ref. 2), is that J(P) is precisely the expected 
score of a particle at P. Thus the right-hand side of Equation 2-34 
consists of: 


f(P) = the direct score at P of J(P) , 

/k(p- I P) J(Pi dP = the score after one or more collision. 


Thus, the function, J(P), represents the importance of particles at P in 
determining < F >. J(P) is commonly called the adjoint function, and this 
term will be used in this paper. The terminology is somewhat ambiguous 
due to the fact that choosing another suitable function other than f(P) for 
Equation 2-34 would have produced a different J(P), and certain authors 
have objected to this terminology (Ref. 3). However, knowledgeable workers 
in the area of radiation transport understand the context where the term 
adjoint is used in conjunction with J(P) and in reference to Equation 2-34. 

2. 3. 2 The Perfect Game 

Since the adjoint function represents the importance of a particle, 
and since * Monte Carlo calculation can be improved by sampling from the 
important parts of the distribution function (i.e. , relative to some effect of 
interest), the question of how the adjoint function can be used to enhance our 
Monte Carlo transport game naturally occurs. To elucidate that question, 
consider the following choices for the parameters in our altered transport 
game: 

(p) . m « -521 (2-43) 

° / K(P'|P) J(P') dP’ + f(P) J(P) 
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Pj(p) 


/ K(P ' | P) J(P') dF' 

/ K{P ' | P) J(P') dP - f(P) 


/K(P |P) J(P') dP' 
J[P) 


§(P) 


Q(P) J(P) s Q(P) J(P ) 

f Q(P') J(P') dP' <3F > 


^(PlP') 


K(P | P ^ J(P) 

/K(P'|P‘) J(P') dP' + f(P') 


K(P | P') J(P) 

ten 

Applying these terms to the transport game illustrated in Figure 2, the 
following parameters can be calculated from the equation whose number 
is given on the left-hand side: 


(2-27) 


w o ( p ) 


= SL?1 

Q(P’) 


<F> 

jrpy 


(2-28) 


Wj(P|P') 


jgPjPl) J(P') . 

K(P|P') W 


As shown in Figure 2, the weight of the transported particle 
is given by 


w 


n 


w 


n-1 


W l (P n! P n-l ) • 


(2-44) 

("-45) 

(2-46) 

( 2 - 4 ?' 

(2-4E' 

(2-49) 
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Substituting Equation 2-48 into Equation 2-49 yields 


w 


n 


w 

W n-1 J[PT~ 
n 


or 




j ( P n ) = 


w n-l 


(2-50) 


Then it follows that 


W n J( p „) = 


(2-5 1) 


From Figure 2, w^ is given by 

W 1 = W o (P l ) { 2 - 52 ) 

Substituting Equations 2-47 and 2-52 into Equation 2-51 yields 

w J(P ) =W (P,) J(P,) = <F> . (2-53) 

n n o 1 1 

Now returning to the transport game, if the particle entered a trapped 
state of the (n + l)th collision, then the score, F, would be 

F . w 0 f(P n )/p 0 (P n ) . (2-54) 

But from Equation 2-43, J(P) can be represented by 

J(P) , JS2L . (2-55) 

P 0 (P) 
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Hence, substituting Equation "-55 into Equation 2-54 yields 


F - w n J(P ) . (2-56) 

But from Equation 2-53, w J(P) ) = <F>, hence 

n n 

F = <F> . (2-57) 

In other words, the exact answer can be obtained from the 
result of the transport game with only one particle. All that is required 
is that the particle be transported until it reaches a trapped state. The 
variance then is obviously zero, and thus, the game played was a perfect 
game. 

However, consider the assumption made to get the values of 
p Q (P), p, (P), Q(P), and K(P|P'), that is, that J(P) is known. It was also 
assumed that <F>was essentially known in Equation2-45, because if J(P) is 
known, then since Q(?) is also known or can be easily calculated, < F > 
can be found by integrating Equation 2-25 by either analytical j r numerical 
technique (i.e. , Monte Carlo is no longer needed). Obviously, if one 
desires to solve for <F> by forward Monte Carlo, then J(P', will prob- 
ably not be known. One possibility is to replace J(P) by some approxi- 
mate importance function, I(P). Then the altered probability distribu- 
tion functions can be found by 


p (P) = — t. (2-58) 

0 j K(P'|P) I(P ') dP ' + f(P) 

Pl (P) = 1 - p Q (?) (2-59) 
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( 2 - 60 ) 


fV 

Q(P) 


UP) Q(P) 

y UP') Q(P') dP' 


K(P j P ') = K(P|P^J(P) 

I(P’) 


( 2 - 61 ) 


Another possibility to compute the altered probability distribu- 
tion functions is to solve for J(P) and ^(P) iteratively, using the values of 
one to define better distribution functions for the other, since if 'A(P) were 
known exactly, a perfect game could also be constructed for the adjoint 
Monte Carlo solution. This procedure is called the iterative forward- adjoint 
Monte Carlo method, and the theoretical validity of this method has been 
demonstrated above. 

2.3.3 Importance and Bias Samplin g 

In the previous subsection, the concept of an importance function 
was introduced, including the use of that importance function to alter a 
phy ically derived distribution function with the new distribution being 
sampled to construct the history of our particles. Altering the distri- 
butions functions, that is, the kernel, K(pj P’), and the first collided 
source, Q(P), lead to much greater efficiency of our calculation for the 
given transport game. This process, known as importance sampling, was 
discussed using the familiar quantity, ^(P), the collision density. For the 
remainder of this discussion of importance sampling, the emergent particle 
density, Y(P). will be the quantity of interest. Although this quantity is 
used less frequently in nuclear engineering, it is much easier to work 
with in Monte Carlo. Thus, the integral emergent particle density equation, 
as given by Equation 2-19, can be written: 
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X(P) = S(P) + f CT(P|P’) X(P') dP' , 


(2-62) 


where 


CT(P | P') - CCeIE’; x) T(x|x’; ¥’) . (2-63) 

Note that the new kernel, CT(P| P’) is not the same as the kernel, K(P|P')> 
defined by Equation 3-20. CT(P|P') requires that the transport kernel be 
applied before the collision kernel, a process which is reversed in K(PjP’). 

A given effect of interest (X), such as energy deposition, biological dose or 
particle flux, can be calculated by: 

X = f g(P)x(P)dP , (2-64) 

where g(P) is the response or payoff function for particles emerging from 
a collision at P in phase space. 

Now consider an arbitrary, but positive, function, I(P), which 
shall be called an importance function. Multiplying Equation 2-62 by I(P)/N, 

where N is the normalization factor of the altered source, given by 
f I(P) S(P) d(P), yields: 


xtP) = §(P) + f x(P') CT(P|P')dP' , 

(2-85) 

where: 


V(P) = x(P)I(P)/N , 

(2-66) 

3(P) = S(P) I(P)/N , 

(2-87) 

CT(P|P') = CT(P| P') I(P)/KP') . 

(2-68) 

Hence, x can be evaluated by: 
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A 


/ g(P) x(P)dP , 


( 2 - 69 ) 


with: 

g(P) = N • g(P) / I(P) . (2-70) 

It does not appear that any improvement has been made by 
calculating X by Equation 2-69 than by Equation 2-64, since only a change 
in notation has occurred. But in the previous section, the use of the adjoint 
as the importance function resulted in a game with zero variance for one 
history. The*. uore, it seems reasonable that a choice of I(P) that approxi- 
mates J{P) should lead to reduced variance. Many studies support this con- 
clusion (see References 1-7). 

An alternate technique which is similar but is based more on 
physical intuition is that of bias sampling. In bias sampling, an altered 
distribution, such as the source distribution, S'(P), or the kernel CT'(P|P') 
is used on the basis of an understanding of the physics of the problem and 
how it relates to the mathematical procedure. In the Monte Carlo game, a 
source particle is selected from S'(P) and then weighted by the ratio of S(P) to 
S'(P). Transport and collision parameters are based on the altered kernel 
CT* (P|P'). Letting X*(P) be the weight density of particles emerging from 
collisions and sources at P, then, 


X '(P) = S'(P) • W (P) + / x'(P') CT(P|P')* W(P|P')dP' 

° J (2-71^ 

where 

W Q (P) = S<P)/S'(P) , (2-72) 

and 

W(PlP') = CT(P| P’i/CT' (P|P') . (2-73) 
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Substitution of Equation 2-72 and 2-73 into 2-71 yields: 

Y'(P) - S(P) +/V(P) CT(P|P*) d? , ( ? -74) 

which proves that \{P) and y'(I , are identical, hence a solution of Equation 
2-71 is also a solution of 2-62. The only restrictions on our altered dis- 
tributions is that all mathematical expressions must be defined over the 
entire phase space, with the indeterminate, ^ , defined as 0. 

The techniques which have been developed ir, this research employ 
aspects of both importance and bias sampling. This is done by using the 
importance function to alter the source energy and angular distribution func- 
tions, the transport kernel, and the collision kernel, then correcting the 
weight of the particles at each step as in bias sampling. For the forward 
Monte Carlo claculation, this importance function is based on both the 
adjoint function and input biasing parameters which allows the use of the 
researchers physical intuition to enhance the variance reduction schemes. 
The equivalent technique is available for the adjoint mode calculation. Thus., 
the advantages of both importance and biasing are included in the techniques 
developed. These techniques are discussed in detail in Sections 3 through 5. 

2. 4 Multigroup Integral Transport Equations 

The computer code from which the test bed for determining 
the validity of the iterative forward-adjoints variance reduction tech- 
niques was derived is MORSE (Ref. 11, 13, 14). MORSE is a multigroup 
Monte Carlo transport code, which means that multigroup cross section 
are used. Because of the strong interaction which occurs between the 
forward and adjoint calculations in the variance reduction techniques, 
an understanding of the multigroup integral equations which are be- 
ing solved by Monte Carlo is essential. In this section, these equations 
are discussed in considerable detail. The beginning point for this dis- 
cussion will be Equation 2-2, the Boltzmann transport Equation. The 
primary emphasis will be placed on those areas where this research 
has deviated from the efforts of other researchers such as Irving in 
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Reference 10, Straker, et.al.,in Reference 6, and Solomito in Reference 
15 and how the quantities used in the importance function are derived. 

For a complete derivation of these equations, Appendix A of Reference 
6 should be consulted, although certain exceptions will be noted in this 
section and Appendix D to that work. 

The time dependence of radiation transport is handled very 
straightforward in Monte Carlo calculations. Each particle is given an 
initial age at the source and the age at some subsequent time is calcu- 
lated by a knowledge oi the particle's speed and distance traveled. 
Therefore, consider the time -independent integro -differential form of 
the Boltzmann transport equation: 

<o . V 0 (x,E, U>) + £ t (x, E) <P (x, E, <3) = S (x,E, Z) 

+ JJ l s (x,E')^> (x,E\i') f(E,5 | E',w') dE’ dw’ (2-75) 

where all terms are the same as defined in Section 2-1, except for 
the assumption of steady state conditions. The multigroup form of 
Equation "-75 is obtained by tegrating over pre-selected energy 
intervals, \E g , where: 


AE 


g 


E - 

g 


E 


R+l’ 


the energy width of the gth group with the highest energy group being 
defined as group 1. Obviously, the sum over all groups, say from 
1 to N, must be identical to the energy from 0 to E^. The resulting 
equation is: 

w.Vrt (x,i) + xf(x)^» (x, «) = S (x, «) 



g ’ JL - - 

X 8 (i)^x,w) 




(u» j u)') del)' , 


(2-76) 
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where. 


(x, li) 



<P (x,E,w) 


dE, 


(2-77) 


r the multigroup (or group) angular flux for group g; 

E 

r K 

S (x,<t>)=J S(x,E,o>) dE, (2-78) 

S — 

E , 

g+1 


= the group source for group g; 




x,E,^)<£ (x, E, o>)dE , 


<p (i,«) 

n 


(2-79) 


energy-averaged total cross section for group g: 


(2 g (x) is defined similarly), 
s 

and 


/ v / 

Vi V» 


f(E, u | E',w')0(x,E,w)dEdE' 


, (2-80) 


/ V *(x f E' f S>) dE' 


E , 

K + 1 


= group g'to group g transfer probability. 


In Equation 2-76, the summation from g' allows both upscattering into higher 
energy groups and downscattering in lower energy groups, although only 
downscattering is important for the problems being considered. 
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The right-hand side of Equation 2-76 is just the expected 
number of particles leaving source or collision events (per unit volume 
solid angle and time) about X and (I whose energies lie in energy r;roup 
g. This is the definition of the group emergent particle density, (x,<n), 
which yields the following identity: 

v r e'- . e-g' 

X (x, o») - S (x,co ) +ZyJ (x)0 _,(x,o/>t fd>U/)dw' (2-811 

B 6 g 6 

2. 4. 1 Forward Multigroup Integral Equations 

Transformation of Equation 2-76 into integral form is 
accomplished by expressing the leakage term in terms of a spatial 
variable, R (as illustrated in Figure 3): 

x' = x - R (2-82) 

where x is a fixed point in space and x' is arbitrary. Taking the 
total derivative of <£g (x,w) with respect to R yields: 


d <p 

<nr g 


(x,^>) 


d#>g dx d<£g dy + d<f> p dz 

dxT dF dy dR dz <3R 


= - rns a d< ^g - cos fl d £ g - COS Y Itjj. . 
dx dy dz 

where a, /I, and y are the angles between dR (which is along o> ) 
and the x, y and z axis, respectively. Thus: 


JL <f> (x,o>) = - w . (x',w) . (2-83) 

dR g g 

Substituting the above equation and (2-81) into Equation (2~76) for 
the point, x' , yields: 

■akV*’") +i t ^ V*^ =x g * 2 " 84 * 
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Introducing the integrating factor , e o' *"t ^ , into (2-84' 

then multiplying by dR and integrating from R=0 to R= oo results 

in: 

co f g _ 

C I,(x-Rw)dR' 

0(x,Z>)=| e" o 1 *g(x -Rc~, Z) dR, (2-85) 

o 

oo r^p . y 
f _ I z (x-R tit) dR\ 

or 0 g (x,o») =1 e g t S g (x-R U) J 

o 


r g , 

+X J Is (x - Hw) <£„,(*- Rw,o>') I )<£’[ dR. (2-86) 

g g > 

Equation 2-36 is called the multigroup integral flux equation or 
the integral flux density equation. Equation 2-85, and hence 
(2-86), can be generalized by performing o :ir integration over 
all spatial points, x', by using the properties of the Dirac delta 
function (see Reference 16) and the differential volume element 
dx' about the point x' : 

dx' = R 2 sin0Wd/dR =R 2 dwdR', (2-87) 

where R - ; x - x'j 


thus: 




(x”)ds 


* g (i', If) 


dx'. 


R< 


( 2 - 88 ) 


where the aneular direction vector from anv arbitrarv x'to the fixed 
point x is defined as a>'(note the change in definition) and the integral 
from x'to x’is along the straight line path, s, which contains the spatial 
points, x? This integral will also be represented by Jx^(x-R'i') dR'. 

0 


From Equations 2-88 and 2-81, a multigroup integral 
emergent particle density eouation can be constructed in a form 
analogous to Equation 2-19. 
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* g (x,w)= s (x,; 


g 


ls(x) 

2f(5) f 


8-g; 


(to | (J) g ((!)'. (J ' 

R^~ 


1) 


■sf(x) 


i 


R 


v#’ -- 


R* w' )dR' 


y .(x.wMdx’dw' 


(2-69) 


2 

where dx' is defined as R du.'"dR (see Figure 4). Tbf group dependent 
kernels are just: 


x) = 


t(x) 


i* 


(2-90) 


and 


T S (k\ x\ to’) = 


8(a< » - 1 ) 

R 2 


if (x) e 


_ jx£(x-R<L' 


)dR 


(2-91) 


However, because of the difficulty in maintaining a consistent 
notation set due to the change of variables of integration for 
different integral equations (such as the flux density and the 
emergent particle density equation) subsequent equations will 
not use the kernel notation unless that meaning is completely 
unambiguous. Thus, Equation 2-89 will be written: 


x e (x, i) = S_(i, u) +T /.Sa. feyjtol 
g 6 if® 


r °° 

. Jlf(x) c“o 


f *V_ - 
_JX| (x-R'^W 


)C^x-Rwi)d«d P , (2-92) 
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with the scattering cross section and denoted by the term, 

(-.-!» * 

It is also possible to write the multigroup collision density 
equation, as shown below. Defining the multigroup collision density by: 


4> (x,o/) = £ g (x) (x,o/ ) > 

& ® 

then from Equation 2-85, and using Figure 4, 

R 


(2-93) 


r co f£ g (x-RV)dR' 

v>(x,w')=J (x) e'o t 


o t 
R 


S (x-Rw')dR 

'£ 


./jf (i) e -A l (i .»;•>«■ I f £ " v®-*"*"* ■ 


(2-94) 


The first term on the right-hand side is the first collision source in the 
g-th energy group and the second term is the group form of the KfPlP*) 
kernel, where P is defined as (x, E*) and P' is (x\ E*) . The similarity 
to equation 2-20 is obvious. The u.oe of the double prime on the angular 
variable is necessary to maintain the geometry notation shown in Figure 
4. However, should the vector entering a collision at x from x' be rede- 
fined ns at' then a>* can be redefined as o> , and the resulting equation is 
then in a more standard form. 

2. 4. 2 The Elfect of Interest 

In Section 2. 3, the solution of Equation 2-25, < F > , could be 
estimated by performing a Monte Carlo transport game based on the 
collision density, ’/'(P). < F > was identified as the effect of interest. 
The effects of interest is some quantity such as the Henderson tissue 
dose at a detector, the total energy deposited in a given volume, or the 
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enei b , flux through a bounded surface. Equation 2-25 shows one way 
in which the effect of interest can be calculated. The response function, 
f(P), is the payoff per unit collision density at P. In multigroup notation, 
the response function shall be denoted by R^(x,o>), which is defined as 
the response function of the effect of interest due to the particle collision 
density at x and in w for group g. 


Due to the simple relationship between the flux and the collision 
density (Equation 2 -4), the payoff function for a unit angular flux at x in 
for group g, 4 > (x,<u), is just: 

D 


R^(x,w) = if (x) Rp(x, J>). 

D D 


( 2 - 95 ' 


The response function is most often given for the unit angular flux, 
since it is usually normaliied to the number of particles per centi- 
meter squared (e.g. rad/ (neutron cm 1. For general use, the response 
function can be a function of position, energy, angular direction and time. 
However, the response function is usually independent of the particle 
direction and age, with the response being defined as zero everywhere 
except at the detector, which could be a point, surface or volume. The 
effect of interest for group g can be expressed as: 


X g = JJr* (x,a>)^> g (x,<L)dw dx, (2-96) 

with the condition placed on X that the total effect of interest, x, 
is the sum of the group effects of interest. 

Group response functions have been defined which pro- 
duce the group effect of interest when integrated over the spatial 
and angular (and time, if desired) variables for both tnc collision 
density and flux. A group response function, R (x,«d), can also 

O 

be defined so that: 
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( 2 - 97 ) 


X 


g 



This is accomplished by requiring that the group effect of interest 
due to the emergent particle density be the same as the effect of 
interest which results from the flux caused by X^(x,o>). Since the 
group flux caused by x^(x,w) at some point a distance of R from 
x along u) is given by: 

r R 

_! if(x+R a))dR' 

0g(x+R <o,cu)= e 0 X (x, to)’ 


( 2 - 98 ) 


C ao , 

R X (k,Z>) X (x,6 >) = J Rg (x+Ra>,a>) <f> (x+ Ra>, oj) dR 


g' ' ' g' ' ' •> * 

0 


f °° J 2t (i+R'o )dR’ ^ _ _ . 

= J e 0 Rg (* + R w, (d ) Xg(x, 61 ' 

0 


( 2 - 99 ) 


from which the definition of R * (x,&») is shown to be: 

o 


Rg(x,a>) = | e‘J 5?(x+ R' cj)dR R £( i+R ^ ^)dR. 


(2-100) 


X - - 

It is not likely that R (x, <•>) will be used directly to calculate X , 
but it will be required in the next section for defining the multi- 
group integral transport equation which is adjoint to Equation 2-92. 
Alternate methods of calculating the total effect of Interest (X or 
< F >) will also be discussed. 
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2. 4. 3 Adjoint Multigroup Integral Equations 

In this section, several forms of the multigroup transport 
equation have been introduced, including the integro-differential 
Boltzmann equation and the integral collision density, emergent par- 
ticle density and flux equations. Likewise, several forms of the ad- 
joint multigroup equation exists, each determining a different adjoint 
function even when the relationship between the forward and adjoint 
"source" term is based on the effect of interest calculation (i.e. A = 

JV (P) f (P)dP =| S(P) J (P)dP Any of these forms can be used 
to calculate the required importance function, but the forms that 
should be used is that adjoint integral equation which is simulated 
in the Monte Carlo random walk procedure in our test bed in the 
same manner as the emergent particle density equation. It is some- 
what surprising that the desired adjoint equation is not the adjoint to 
the emergent particle density equation. For that reason, this sub- 
section contains a discussion of different adjoint equations and the 
relationship between these equations which determine the importance 
function. Hue to the adequate descriptions of the derivation of most of 
these adjoint equations in References 10 and 11, a full derivation of all 
equations will be omitted. However, errors, other than typographical, 
in Reference 11 will be discussed and corrected in the following develop- 
ment or in Appendix D. 

One method of defining an adjoint equation is to derive the 
adjoint to the inter gro -differential Boltzmann equation (2-75) by the re- 
quirement that: 

J <M P) L 0 (P)dP = J 0 (P)lV*(P) dP, (2-101) 

* 

where is usually called the adjoint flux and L is the Boltzmann 
transport operator, defined by: 


2-36 



L<*>(P)= - £ -V 9 (i, E) - X t (x, E)*(x, E) 

+ J J l s (x,E ) f (E,£>| E\«')*(i',E')dE*du. 

The resulting adjoint operator, L* , differs from the forward 
operator in two ways: 

1) The leakage terms have opposite signs, 

2) The transfer probability in the inscattering terms 
have reversed the initial and final direction (to) 
ai’d energy (E) variables. 

L/V*(P) = w . V0*(x, E) - I t (x, E) 0*(x,£) 

+ J J Xgfx.E*)! (E'.i'lE.iWx.EO dE'dw. 

Defining S* (P) as the adjoint source term, the multigroup form of 
the adjoint equation, 

L*0*(P) - S* (P) = O, 

can be determined as was done for Equation 2-75. If the adjoint 
flux, 9 *(P), is used to define the multigroup cross sections and 
transfer probability, then the cross sections for the adjoint Monte 
Carlo calculations will probably differ from these for the forward 
calculation. However, if the same multigroup cross sections are 
used and only the transfer probability matrix is transposed, then 
the multigroup form of Equation 2-104 is identical to the adjoint 
form of Equation 2-76: 

- w . \ £g (x,a>) + (x)4>*(x,w) = Sp(x,«) 

+ X [ if B (i;&»|w)* (*,«') 

r J g 

Remember that the choice of S* (x, Z> ) will determine the solution 

o f (x, to ) . 


( 2 - 102 ) 


(2-103) 


(2-104) 


(2-105) 
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Using the same procedure as applied to Equation 2-76, 
but with x' = x+ Rw and an integrating factor of: 

" Jlf(x+R'w )dR' , 
e 0 

the following integral equation is derived: 

* f 00 I ^ g 

V*’“ )= J o s‘(it R ;,oidR 

0 

°D f R e 

+ j I? (x + RL) e~J + R^)dR' 


*?< 


r g* _ g, - - * 

I i s (X + Ro> ;ai [ o ) • 0 

g' ^ 2 ^ (x *■ Rcj) 


(x+R a*, a> ) dw‘ dR. 


(2-:o:>) 


The above equation is adjoint to Equation 2-92, the emergent particle 

density equation, since the second term on the right hand side is just 

the adjoint (or transpose in th s case) of the corresponding term above. 

The above equation was derived from the integro -differential equation 

and is not the adjoint of the integral flux equation (2-86). Because of 

the ambiguity of terminology introduced by 0* (i,L ) , Equation 2-106 

* _ _ 8 
will be denoted by X^(x, « )• That is: 


X* (x,w) = ** (x,i) . (2-107) 

g g 

At this time, the source or first term of Equation 2-106 is unspeci- 

* 

fied (which means that X (x >t L) is unspecified). As discussed in 
Section 2. 3. 1. the choice oi the source term should be the one 
which allows the effect of interest to be calculated by either the 
forward or adjoint function. This means that since: 


A = 



X (x, u>) R*(x, w)d<jdx 

*> D 
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(2-108) 



f J Y (x,Z* ) S (x, l , ) doidx , 

J * U P 


then the adjoint source term is: 

f 00 f ^ 6 , 

R^(x,w) = J e”J (x+R‘w)dR' s*(x+R<L,Zj )dR. (2-109) 

0 

Equation 2-106, which is called the integral point value equation, 
can be written: 

r R 

X*(i,L) =R v (x,w)+ t lf(x+Ri)e'J if (x+R'a>) dR' 
g g J t 0 



V ^ r> “ 

Ig (x+Rcj 


V ? (x- Ra>) 


— - y g .(i+R«,i')d«dR. 


( 2 - 110 ) 


Comparison of Equation 2-i09 with Equation 2-100 provides the 
identification of the adjoint source term, S* (x,a»), which was 
first used in multi-group integro-differential adjoint equation (2-105): 

S = rJ(*,w), (2-111) 

o o 

the angular flux response function. It is this function that should 
be input as the adjoint source into computer programs such as 
ANISN (Ref. 17) and DOT (Ref. 19) which solve the multigroup ad- 
joint integro-differential equations using finite difference techniques. 

Equation 2-110 represents one form of the multigroup 
adjoint integral transport equation. While this form is solvable 
by Monte Carlo methods, it is not solvable by the same techniques 
used in the test bed to solve the forward emergent particle density 
equation (2-92). This is obvious because the first step alter the 
source particle selection is the collision kernel step, whereas that 
step in simulation of x g (x, w ) is the transport step. Therefore, 
consider the adjoint form of the collision density equation (2-94): 
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U> (x, hi) - R'*' (x 
B B 


.00 




g — g 

>-S ( X ’ (U I <*j ) 


R / 


If S (x) 


r re' , * 

vK' ( - + R - 0e -J (x + F cj ) dR' ^ p (x ^ Rw'.u) dR dd-\ (2-1 


12 ) 


This equation is the multigroup form of the adjoim equation for 
J(P) derived in Section 2. 3. 1. It is also called the value, im- 
portance, and integral event-value equation. The designation in 
this paper will be the value equation, because the equation can be 

derived by defining «J> (x,o> ) to be the value of an event or collision 

B 

at x to the effect of interest for a particle in group g entering a 

collision with direction cj . is the immediate payoff and the 

B 

second term represents all subsequents contributions. It can also 

be shown that X * and are related by: 

S B R 

0 


^*(x + Roj, »)dR. (2-113) 

* 

Thus if is the value for a particle upon entering a collision, 

* * 

then the above equation indicates that is the value for a particle 
leaving a collision. That this is the case is also indicated by the 

y 

fact that the source term in Equation 2-110 is R , the response 

B 

function due to the emergent particle density. 

Inspection of value equation indicates that it is in the 
same form as the emergent particle density equation (Xg(x,o*)) ( 
which means that it may be simulated by the same Monte Carlo 
methods. Hovever, implementation of the simulation for the value 
equation soon experiences difficulties. After selection of an "ad- 
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joint particle" for x|w',and g' from P^(x,w), the adjoint particle 
is not only traveling in a direction opposite to the direction <L (i. e. 
from x + R wto x), but the transport kernel is not normalized, since 
(with x + Roi = x ): 


i-f (x+ Ro7) 



R «) dR' 


= *1 (iO 



and our object is to pick a point x = x* - Rtc'. This obstacle can 

be overcome by choosing from a properly normalized kernel and 

correcting the weight of the adjoint particle by if (x\u)/ if a , Cl) ) . 

However, this calculation is not required in the simulation of V (x, i>), 

S 

so it would be better to avoid the correction in the adjoint equation 
simulation. In addition, implementation of the collision kernel reveals 
that it is also unnormalized, and that the corrections required are 
even more extensive and time-consuming. Our problem with the 
transport kernel can be circumvented by defining a new function which 
is the product cf </^ and if, since the introduction of if (x) 
produces a properly normalized transport kernel. 


The new function will be denoted by G (x,w) and will be 

defined by: 

G (x,uj) =lf(x) ^ (x,o») 
g 1 g 

Substituting the equation for (x,u>) into the above equation yields: 

g'— g 

G g (i,i) = sf(i)R g U (x,i) + Xj if(i) 


( 2 - 114 ) 




00 f^g 

f e -J i t (x+Rw)dR kf(i + Rw ) 0 * (i+Rw'.i’jdRdi’, 

J 0 o 


(ri * j m 
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which can be written as* 


G,(i. i 


= K*<x.£i> + X/ 

g’ 


if " g (x;a'’| <L) 


R 

00 r 

’[ lf'(x)e"/ lf'(x t- R'cc') dR’ G lx +Ri’,cii' dR uu:' 

' 7 l ^ L cr T 

O & 

° _ 

Equation 2-95 was used to define R^(x, u>) and the unity ratio, 

cr* ~ . <t’ _ ts 

i. j (x)/ (x), was distributed in the second term to normalize the 

transport kernel. The above equation also allows the adjoint particles 
to travel in a direction opposite to their velocity vector. To overcome 
this somewhat confusing convention, define the adjunction, an "ad- 
joint particle” that obeys Equation 2-116, except that the directions 
have all been reversed. This allows the adjuncton to travel in the 
same sense as the velocity vectors, but also requires special attention 
in the interpretation of the results of this adjoint calculation, since all 
results will be in the opposite sense to the value calculated by Equation 
2-112. The new equation defined is: 


G (x, w) 
6 



if (x) 


CO R 

/lf'(x) e f sf ,(i - Gg/x - da?' 

o (2-117) 

Equation 2-117 will be called the integral emergent adjuncton 
density equation, which is in agreement with Appendix A of Reference 11. 
However, the approach taken here is different from that taken in Reference 
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11, an approach which is patterned after Irving (Ref. 10). The deriva- 
tion of the emergent adjuncton density equation contained in Reference il 
has several errors, mostly in terminology (e. g. ,00 distinction is made 
between G and G ). Appendix D contains a corrected derivation of the 
multigroup emergent adjuncton density equation which is consistent with 
the definition of (x,7>) in Reference 11. The result ir. Reference 11 is 
similar to Equation 2- : 17, but it requires considerably more effort than is 
required by the derivation of G (x,u?) based on its relationship to the value 
equation. Comparison of the emergent adjunction density equation (2-117) and 
the emergent parcicle density equation (2-92) shows that the Monte Carlo 
method used to simulate the emergent particle density equation can be 
used to simulate the emergent adjuncton density equation. Only the source 
functions, S (x,a») and R^ (x,o>), and the scattering kernels are different. 

A description of this simulation will be given in the next section. 


2. 4. 4. Effect of Interest and Importance Function Estimation 

Now that the integral transport equation, both forward and 
adjoint, which will be simulated have been chosen, two other questions 
must be considered: 

1) Hov; will the effect of interest be calculated? 

2) What functions will be used in determining the impor- 
tance function? 

Most of the necessary information for considering these two questions 
have already been developed in this section, and only need to be sum- 
merized at this time. 


Consider first the relationship which was required be- 
tween the forward integral equation and Its adjoint, that is, the 
adjoint function must not only meet the mathematical requirements 
but also that the adjoint source term must be the response function 
used to generate the effect of interest in the forward mode. This 
uniquely related the following sets of equations: 

. collision density (4» ) - value (»|£) 

. emergent particle density (x ) - point -value (V ) or ad- 
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joint flux ( 4> ) 

. flux density (0 ) - emergent adjuncton density (G or 
E E 

G ) 

£ 

Since the most natural method of calculating the effect 
of interest, A, is: 


* III R^(x,aO <f> (x,oj) doidx, 


and the source term for <P (x,<L)from Equation 2-86 is: 

s 


r oo r g 

J e .jv t (x-R'd>')dR / g (x - R o>, oj)dR, 


then the effect of interest can also be expressed as: 


■ZfJIfV 
£ 0 


00 - I I,(x- R’Z>)dR 


S (x -RtI>,tu)dR 
g 


* G^(x, u>) dwdx . 

Obviously, the above expression is not a very desirable method 
to calculate \ . But since the emergent adjuncton density equation 
is being simulated, and the flux density is not, won't it be necessary 
to use Equation 2-119, instead of 2-118? The answer is fortunately 
no, since in the forward Monte Carlo random walk, each of the quan- 
tities, 4* , X and 0 , can be estimated if desired. Likewise, al - 
g £ £ 

though the emergent adjuncton density, G’ , is being simulated, it 

* * ^ 

is also possible to estimate X and 4> . While this procedure will 

g g 

be discussed in the next section, equations such as 2-85, which re- 
lates <£ to X , and 2-114, relating G and provide an insight in- 
g g g g 

to how this is accomplished. Therefore, there is no restriction as 
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to what quantities, either forward or adjoint, can be used in the 
effect of interest extimation. 


The forward estimation of X will probably be with the an- 
gular flux, yWx,d), in the Monte Carlo approximation of Equation 
2-118, because the response function is the one most often 

use«' *r radiation transport calculations. For the adjoint estimation 
of A, the adjoint response function most likely to be available is 

the source function, S (x, d). This means that the adjoint flux, 0 (x, d) 
* & __ 5 

or X (x, co ), will be determined in the G .(x, a’) simulation and 

o 6' 

the contribution to effect of interest calculated. This estimation is 
based on the fact that: 


■ dd I 1 S (x,d)X (x, d ) dx d d. 
g J 1 g g 


( 2 - 120 ) 


It should be recognised that while the group effect of in- 
terest, A or A* , can be defined for each group In both the forward 

n P 

and adjoint calculation, these values are not necessarily equivalent, 
contrary to the many equations in Reference 11 which state that they 
are. This fact r*b» demonstrated by considering a two group 
problem, vhti jwinw are emitted in the first or fast group 
and where the defector haa a non -zero response only for neutrons in 
the second or thermal group. Assuming a infinite, non-absorbing sys- 
tem. then the forward and adjoint flux will be non -zero for both groups, 
but the sdjoint uource (or forward response) is zero of the fast group. 

If the thermal group effect of interest is evaluated by the forward flux, 
then: 


= J J R^(x,w)d> T (x, Zj) dddx 


> 0 , 


0 

since both R T and 0^, are positive definite values. However, if the 
thermal group effect of interest is evaluated by the adjoint flux, then 
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* 




- o. 


since the thermal source, S T , is identically zero. Therefore, A T 
is not equal to A* , and the forward and adjoint group effect of in- 
terest cannot be compared directly. 


Selection of the functions to be used for determining the 
importance is a much more different problem than that of the effect 
of interest calculation. However, it was shown in Section 2. 3. 2 that 
the adjoint or value function, J(P), could produce a perfect game for 
the specified transport game which simulated the collision density 
equation. This result was investigated even more thoroughly and 
extended to other transport game and biasing schemes in Reference 
3. The authors concluded that the value function is always a good 
choice for importance function biasing of the collision density equa- 
tion. Therefore, it was decided on the basis of the conclusions 
above and because the effect of interest estimation required the 
same functions that the following functions would be used in the 
importance function: 

♦ _ _ 

Forward; The adjoint flux, X (x,cj), for biasing the emer 
gent particle density equation simulation. 

Adjoint: The angular flux. 0 g (x,tL), for biasing the emer- 
gent adjunction density equation simulation. 

This procedure is illustrated by the diagram below: 


Mode: 

Forward: 

Adjoint: 


Equation Importance Function 

Simulp.t' ’ Estimated 


* g (x,w). 


G r (x,w) - 


<? g (x,w) 


-X*(x,Z>) 


The arrow from the "Importance Function Estimated" column to the 
"Equation Simulated" column indicates the equation in which the im- 
portance function is used. The exaci form of the importance functions 
will be discussed in Section 4. 
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3. THE IFAM TEST BED 


The principal objective of this research is to develop and evaluate 
variance reduction techniques which utilize the iterative forward- 
adjoint Monte Carlo method. A necessary tool for accomplishing this 
objective is a computer program which utilizes these techniques while 
performing the Monte Carlo simulation of the ii egral transport equations. 
No such computer program existed, so one had to be developed, coded 
and validated as part of this research effort. Because this program 
must be capable of accepting several different importance and biasing 
functions, it was decided that the proper approach would be to develop 
a computer software test bed. The test bed must handle, for both 
forward and adjoint modes, all input and output of data, cross section 
manipulations, tracking of particles in complex geometries, and 
provide the framework for the Monte Carlo random walk, importance 
function estimation, and effect of interest calculations. A test bed 
is distinguished from a computer program or code by the fact that 
the test bed requires the addition of subprograms or portions of a 
subprogram which perform specific functions (in this case, the 
particular variance reduction schemes) before it is a complete and 
executable program. 

A test bed is written in such a manner as to minimize the impact 
of adding the coding required to perform these functions. This is neces- 
sary since a major consideration in the evaluation of the different 
algorithms is computation time, and the time required for computing 
activities independent of the algorithms should be constant. Also, the 
testing of a particular algorithm should require only the validation of 
the coding for that algorithm, and not the validation of the entire 
software package. These requirements were implemented in the 
development of IFAM (Iterative Forward -Adjoint MORSE), the name 
of the test bed utilized in the evaluation of the variance reduction technique. 
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This section contains a description of the IFAM test bed. This 
description is necessary for an understanding of the 'evaluation of the 
different techniques which will be discussed in Section 4. In many ways, 
the test bed represents to the analyst what equipment and instrumentation 
represent to the experimenter. Therefore, the IFAM description has 
been placed in the main body of this paper, beginning with a brief 
description of the computer code MORSE (Ref. 11), which forms the 
basis of IFAM. The rationale for selecting MORSE instead of other 
codes as the IFAM basis is also discussed. The simulation of the 
integral forward (X (x, cc)) and adjoint (G (x, w)) transport equations 
are also described, followed by a discussion of the test bed logic or 
flow and a short description of IFAM-peeuliar subroutines (e. g. , FAIF, 
which is used to calculate the importance and biasing functions). 

3. 1 The MORSE Code 

Since there was no computer program which performed the iterative 
forward -adjoint Monte Carlo radiation transport calculation, development 
of the test bed wa& essential. One appraoeh was to develop and code 
the test bed from "scratch”, except for the cannibalization of certain 
parts of other transport codes. However, in order to include the 
capability to solve a wide range of transport problems in complex 
geometries, this approach would require many man-years of effort. 

A much faster approach is to utilize one of the already developed 
Monte Carlo transport codes and make appropriate changes so that 
the iterative forward-adjoint calculations can be performed. Those 
codes which were considered as the basis for the forward and/or 
adjoint Monte Carlo calculations include MORSE (Ref. 11), CAVEAT 
(Ref. 19), 06R (Ref. 20), 06R-D (Ref. 21), SAM-CE (Ref. 22), 

ANTE -IT (Ref. 23), and GADJET (Ref. 24). The criteria on which 
the code selections were based included speed, accuracy, representa- 
tion of the physical phenomenology, geometry handling schemes, 
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input and output compatability between the forward and adjoint modes, 
ability to perform both forward and adjoint calculation, and applicability 
to the iterative forward -adjoint method. Also, the author’s familiarity 
with the code and the effort requ red to modify into the test bed were 
considered. Comparison of the above codes soon revealed that MORSE 
was far superior as the basis for the lest bed. 

MORSE is a multigroup code which handles neutron, gamma ray, 
and coupled neutron-gamma ray problems for both the forward and 
adjoint mode. The SAM-CE code system can perform forward neutron 
and both for\ , 'ard and adjoint gamma ray problems. However, SAM-CE 
requires three separate programs to accomplish these calculations and 
it does not have a provision for neutron adjoint calculations. CAVEAT 
is also a code system which handles the forward calculation of neutron, 
gamma ray, and coupled problems, but performs no adjoint calculation. 
Thus both CAVEAT and SAM-CE would require additional code develop- 
ment to completely handle the adjoint calculation. The other code 
candidates are restricted to either neutron (06R, 06R-D, ANTE -II) 
or gamma ray (GADJET) problems and only in one mode (forward 
or adjoint). 

Therefore, the MORSE code definitely is better in respect to the 
ability to perform both forward and adjoint calculation and applicability 
to the iterative method. The criteria of speed and accuracy are very 
difficult to evaluate, but it is generally accepted that MORSE is faster 
than the other candidate because of its use of multigroup cross section, 
although comparisons of accuracy (Ref. 25) have been inconclusive. 

Only MORSE is a multigroup code , the other candidates use point cross 
sections. This means that the physics of the collision process is contained 
in the multigroup cross sections in MORSE, but are handled explicitly 
in the other codes. However, any disadvantage which MORSE may acquire 
in the representation of the physical phenomenology is balanced by the 


3-3 



advantage of the multigroup code in the ease of generating the adjoint 
cross sections. 

Of all the geometry schemes studies, the co :ib. natorial geometry 
package, developed for the SAM-CE codes, is best. Since there exists 
a version of MORSE which contains this geometry scheme (Ref. 13), 
then only 06R, 06R-D, and CAVEAT are at a disadvantage to the 
other code candidates. MORSE also shows a very definite advantage 
over the other candidates in input and output compatability, since 
only a few input quantities need to be changed to go from a forward 
to an adjoint calculation. Only for one criterion does MORSE take 
a definite second place, and that is in the author's familiarity with 
the code. At the time the selection was performed, the author bad 
a detailed understanding of all aspects of CAVEAT, but only the 
combinatorial geometry package in MORSE. 

MORSE was also considered to require less effort to modify into 
the test bed. The obvious choice for the test bed basis was MORSE, 
and the test bed name, IFAM or Uerative Forward-Adjoint MORSE, 
reflects that choice. The major modifications required to execute 
this transformation were: 

1) New logic so that both forward and adjoint calculation can be 
performed in the same run; 

2) New logic for inputting both forward and adjoint data and 
the temporary storage of the data for the mode not being 
calculated; 

3) Estimation and storage of the quantity to be used in the 
opposite mode importance function; 

4) Storage and application of the importance function for 
biasing the distribution functions; 

5) Calculation o! the new importance and biasing functions 
between the mode calculations; 
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6) Setting up an overlay (or segmentation) map to increase the 

amount of core storage available for the importance function. 
Other modifications for convenience (e. g. storing the cross section 
input for the forward mode so they can be read from a data file 
instead of the input unit during the adjoint input operations) were 
also implemented. 

In order to understand IFAM, a familiarity with the MORSE code 
is required. Many of the features and attributes of MORSE have been 
discussed above. Other characteristics that are needed to understand 
IFAM are given below. Complete details on MORSE and its analysis 
package. SAMBO, are contained in Reference 11, 13, and 14. In this 
summation, the term MORSE will include the basic code package 
(Ref. 11), revisions to MORSE such as those to the geometry routines 
(Ref. 13), and the analysis package (Ref. 14). 

As discussed above, MORSE is a multipurpose Monte Carlo trans- 
port code for both neutrons and gamma rays. Both forward and adjoint 
problems can be solved with little change in input data since multigroup 
cross sections are used. All large data arrays have been included in 
the variable dimensioning feature to optimize the use of core storage. 
Complex three-dimensional, time dependents problems can be solved. 

In addition, several types of importance sampling options are available, 
including: 

• Source energy biasing, 

• Path length stretching, 

• Downscatter energy biasing, 

• Russian roulette and splitting. 

(These options have been maintained or upgraded with data from the 
importance function in IFAM. ) 
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By employing muMgroup cross sections for each different material 
in the problem geometry, the physics of the interaction or collision 
process is contained implicitly in the cross sections. The type of 
particle, neutron or gamma ray, transported affects only the cross 
sections, so that the logic is the same for both, except for certain 
"editorial” exceptions. The multigroup cross sections account for 
anisotropic scattering by requiring that each group-to-group transfer 
has an associated angular distribution which is calculated from the 
input Legendre coefficients using a generalized Gaussian quadrature. 
Isotropic scattering is handled by randomly selecting the outgoing 
angles or direction cosines. The crnss section module processes 
the input data so as to produce a normalized collision kernel which 
is multiplied by the nonabsorption probability for the forward mode 
and by the analogous but non-physical weight factor for the adjoint 
mode. (This weight factor will be defined later. ) 

The geometry module in the MORSE version used as the IFAM 
basis is the Combinatorial Geometry package (Ref. 13). This package 
describes general three dimensional material configurations or regions 
by performing the union, differences and intersections of simple bodies 
such as spheres, boxes, cylinders, wedges, cones, and ellipsoids, 
and a more complex arbitrary convex polyhedron of four, five, or six 
sides. Experience has shown that almost any configuration can be 
represented in this manner. 

Scoring options include those available in the SAMBO analysis 
package (Ref. 14) plus built-in collision density and track length fluence 
estimators (Ref. 13). SAMBO allows an arbitrary number of detectors, 
energy-dependent response functions, energy bins, angle bins, and 
time bins. The scoring can be divided into: 

• Uncollided and total response, 

• Fluence versus energy and detector. 
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• Time dencndent response, 

• Fluence verses time, energy, and detector, 

• Fluence verb s angle, energy, and detector. 

The standard deviation i.s calculated for each of the above quantities. 
Estimators for point, surface, and volume detectors are available. 

The technique of estimating the above quantities by simulating the for- 
ward or adjoint integral transport equation is explained in the next 
section. 

3. 2 Monte Carlo Simulation of Integral Transport Equations 

Both the forward and the adjoint transport equations are simu- 
lated in the IFAM test bed. The same coding is used for both of these 
simulations, with changes only in input data (e.g. , the source term and 
the cross section scattering matrix) and the identification of the quanti- 
ties calculated. The Monte Carlo simulation of the integral equations 
is discussed below in a very simplified manner. Details of how the 
specific phase space components are selected are contained in Chapter 

4. This discussion emphasizes the mathematical and physical aspects 
of the simulation. Also, explanation of the effect of importance sampl- 
ing on the probability distribution functions will be delayed until 
Chapter 4. 

The simulation of the forward integral transport equation (2-92) 
and the adjoint equation (2-117) will be considered simultaneously. 

Thus it is oossible to compare both the similarities and differences in 
each of the three simulation steps. Table I shows the major elements 
of the simulation. For each of the three steps, the table defines integral 
equation terms which are simulated, radiation quantities which are esti- 
mated, and the expression for the value of the estimate. The treatment 
of both discrete and continuous probability distributions will be 
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Table I. Monte Carlo Simulation Outline 


0 ) 

13 

E 

■r“ « 

vj 

W 


CD 

3 

> 


TJ 

Ol 

-4-1 

OS 

£ 

-M 

W 

w 

U) 

CD 


<U 

-4-4 

J2 

3 


-4-4 

CO 


0> 

s 

s 


' 3 
T3 

IX 

•c 

'*3 

iX 


CO 


u. 


Wc 


bC 


i 3 
iX 


> 


• 3 

iX 

be 


® 

O 

o 

m 


T 3 

<3 

* 

£ 


i 3 
T3 

1 X 
T 3 

• 3 

Q< 

S« be 
■ OS 


£ 

be 

s 

£ 

d 

u 

•H 

•*-1 

3 

(0 


b£ 4_j 

^4 


l 

01 

f—4 

u 

rt 

tc 




W 


be 


M X 
£ £ 


he — w 
W 

• 

oo a 

r* £ 


i 3 

i’x 


i 3 


•a 


be 


• 3,3 

IX 
£ 


• 3 

i 

IX 


I 3 

i 

;*x 


be 

9^ 


* be * be 
x x 


13 

cc 

i 

"be 4_ 
M 


! 3 
a: 

i 

l* 

"be *4 

14 


CO 


'3 
• >< 

be 


i»< 

"be *4 
W 


i 3 

• 3 


£ 

o 

•o 

< 


■ X 

1 

I*x 


u 

a 

a 

X 


T3 

el 

l 

£ 


n 

—■+ 

o 

"O 

*< 


•3 

t 

'*3 

be 

t 

be 


e 

O 

•H 

OB 

a 


£ 


m 


•3 


tx 


I U 


<>< 

be -4-4 

i-i 


'3 

'*3 


E 

• 3 

• 3 



■ X 

IK 

to 


- 

, 


— 

be 

£ 

u 

|-L 

IX 

I^X 

IX 

be 

1 

i 

* 

0 ) 

~b 0 

«i*bo 

be — 

"be 4-4 

be on 

'be a 

H 


• K 

W 

(4 

(4 

14 


3 , 

•o 

< 


3-8 



considered achievable in the sense of the Lebesque - Stieltjes integral. 

That is, the probability of the events between Xj and x^ is just: 

*2 

Pr [Xj < XS Xj] = /dF(X) (3-1) 

X 1 

where X is a random variable with values x, and F(x) = Pr [ X <, x] is the 
cumulative distribution function. Assuming the existence of a set of 
random numbers uniformily distributed on the interval from 0 to 1, then 
a specific value of x can be selected from f\x) by the following relationship: 

x 

P = J dF(X) (3-2^ 

- 00 

where P has been randomly chosen from the uniform set of random number. 
When the above relationship is solved repeatedly for x with p's taken from 
the uniform random number set, then the x's will be distributed as F(X). 
This principle is the basis for the Monte Carlo technique, and it is in 
this sense that the expression "select ..." is used below. 

3. 2. 1 Source Selection 

The source term in Equation 2-92 is expressed as S (x, to). 

8 

This term is completely general, and its exact form depends on the physical 
problem being solved. If our source is an isotropic, mono -energetic point 
source, then the source selection requires that only the energy group and 
spatial coordinates be specified in the input. The angular coordinate, 
usually expressed by the direction cosines, can be selected rather simply 
by assuming an uniform distribution over the unit sphere. Details of this 
selection will not be necessary for this discussion, since it will be assumed 
that the selection of the energy group, position and direction is possible 
from a properly defined source term. Thus, select the source parameters 
(g\ x’, «') from W' 1 - S , (x', w’), where: 

8 D 
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(3-3) 


W S '£,// s d “' di ' • 

O 

The term, W~\ normalizes our source term so that the selection is 

performed on a properly normalized distribution function. W is also 

s 

the statistical weight assigned to each particle whose parameters are 
selected from the source distribution function. 

Since the particle selected above is emerging from the source, 
then the weight of the particles, W is an estimate of the emergent 

5 

particle density, V , (x’ , uj'). At this time in the simulation, W should 
be used to "score" the emergent particle density, if desired. This is 
usually done by keeping a running sum of the statistical weights in finite 
elements of phase space. 

— S-* — — 

The adjoint source term, given in Equation 2-117, is R (x, oj). 
— _ © 

This function of g, x, and os must also be normalized so that these adjoint 

source or adjuncton parameters can be selected. After the adjuncton 

parameters (g*, x’, «’) have been selected, the adjuncton is assigned a 

statistical weight of: 

w * =£// R^, ( x', i')dS' dx'. (3-4) 

As in the forward source selection, the emergent adjuncton density, 

G , (x, go’), can be scored by keeping a running sum of W and other 

adjunction weights as discussed below. However, neither G (x’,w’) 

or X , (x’, Z') are usually scored since they do not yield estimates of 
(5 

the effect of interest when multiplied by the normal response functions 
(e. g. , see Equations 2-97 and 2-100 or 2-119). 

3.2.2 Transport Step 

Inspection of the integral term Equations 2-92 and 2-117 shows 
that the next step in the simulation of the emergent particle density or 
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emergent adjuncton density is the transport of our "particle” from x' to x 
along the u>' direction. The next colli sion site, x, is selected from the 
probability distribution function: 

r R c - 

, - Lf (x - R' ui ' ) dR' (3-5) 

Ef (x ) e J Q 1 

It is possible that the new collision site is outside the region ol 

interest or has escaped into a non-reentrant void. In this case the history 

of this particle (or adjuncton for the adjoint mode) is terminated and a new 

source particle is selected. If x is within the region of interest, then this 

transport step can be used to estimate the collision density, </> , (x, oo’), 

K 

and hence the flux, 9 , (x, oo), in the forward mode. The estimate for 

the collision density is just the statistical weight, W , while the estimate 

t s 

of the flux is W /if (x). These estimates are also made for finite elements 
8 t 

of phase space, with the spatial bin being that geometric region which 
contains x. 

It is also possible to get an estimate of the fluence in those 
geometric regions through which the particle passes in traveling from x' to 
x. This estimate is based on the fact that the fluence can be defined as 
the path length per unit volume. Therefore, the contributions to the estimate 
of the average fluence in a given geometric region when a particle travels 
a distance S through the region is just the product of the statistical weight W 
times the distance, and divided by the region volume, V (i.e. , S • W /V). 

A running sum is usually kept of the product of S and W for each energy 
group and region during the random walk calculations. This sum is 
normalized by dividing by the region volume and the total number of source 
particles at the conclusion of the calculation, producing an est..nate of the 
fluence. Flux estimates can also be calculated for time dependent problems 
by dividing the time range into time bins, recording the running sum estimate 
for each time bin as well as energy and region, then including the time 
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difference in the normalization term. For time dependent problems, the 
flux estimate calculation is identical to tho fiuence calculation, except the 
source term is in units of per second. 

Similarly, the transport step of the adjoint simulation produces 

estimates of the adjoint flux. The statistic 1 weight at each collision point 

can be divided by the total cross section to produce an estimate of the 
* _ 

adjoint flux, (x, - to). Track length estimates of the adjoint flux 
during the adjoint simulation are calculated in the same manner and use 
the same computer coding as the 1 'orward flux. Running sums of these 
quantities are calculated by 1 7 A’! for each energy group, region and 
angular bin. 

3. 2. 3 Collision Step 

After the particle or adjuncton has been transported from either 
a source or a previous collision site to a new collision site (e. g. , from 
x’ to x), the physics of interaction was a nucleus must be simulated. This 
process is usually assumed to occur instantaneously and at the selected 
collision site. Mathematically, the collision step is a simulation of the 
scattering or collision kernel, which results in the selection of a new 
energy or energy group (g* to g) and a new direction (w’ to w). In order 
to make this selection of a new energy group and direction, the collision 
kernel must be properly normalized. This normalization for the forward 
mode can be performed in two ways: 

1) Terminate (or kill) the particle or adjuncton with an 
expectation proportional to the absorption probability. 

2) Adjust the statistical weight W by multiplying by the 

s 

non-absorption probability. 

The second method is the one most commonly used, and is the method 
used in IFAM. 
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This method has the effect of dividing the collision term into 

two parts: 

I g (x) £ g ** g (x . clj ' u’ ) 

s * ’ 

r g ’ (x) r|'(x) 

which is the same as t.e term shown in Table 1. However, in this r,rra, 
the first ratio represents the probability that the collision event will not 
‘■erminate the particle. Reducing the statistical weight by multiplying by 
this ratio has the effects of "killing" that fraction of the particle that would 
"on the average" terminate the collision site. Of course, a part of a particle 
cannot be terminated in the physical process, but this mathematical 
technique is perfectly correct. It then allows the selection of the new energy 
group and direction from a normalized function, as shown in the second 
part of the above term. 

For the adjoint simulation, an equivalent mathematical procedure 
is used to simulate the adjoint collision term. However, the adjustment 
of the statistical weight by 

[E / Eg’ " w'| w M )d^”] /l t g '(J) 

g” 

has no physical analog such as the non-absorption probability. In fact, 
the above term can have a value greater than one, unlike the non -absorption 
probability. Use of the above term does allow simulation of the adjoint 
integral equation for (x, tu) in a manner identical to the forward into ral 
equation for x (x. o>). It also allows the transport step to be performed 
with a normalized kernel, as discussed in Chapter 2. Only the cross 
sections, which are usually input for each problem, differ in the simulation 
of the particle (neutror and gamma ray) and adjuncton transport. 

After the adjustment of the statistical weight to allow normali- 
zation of the collision kernel, selection of the new energy group and 
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direction is performed from the following term: 

j-g - fi (x ; go' | w)/j^J r g ’ g (*'• u.'’ | a.'") da.'”J 

The general procedure for this selection is to determine the downscattering 
probabilities into each group g from the current group g'. These prob- 
abilities are precomputed and stored by the IFAM test bed. Then the 
outgoing angle is selected from the discrete angle probability matrix 
which corresponds to the selected downscatter energy. If the new energy 
group is outside the range of interest (e. g. , below the energy at which 
the response function for the effort of interest is zero), then the particle 
or adjuncton history is terminated and a new source is selected. If the 
history is not terminated, the transport step is performed with the new 
parameters (g, x, uj). 

The estimation of quantities of interest at the conclusion of the 

collision step is very similar to that of the source selection step. The 

outgoing statistical weight is an estimate of the emergent particle density 

Xg (x, ui) for the forward mode and the emergent adjuncton density 

G (5, u>) for the adjoint mode. Of course, the statistical weights can 
S 

change from step to step, especially at the weight adjustment in the 
collision step. For a problem in which importance sampling techniques 
are utilized, weight adjustments will occur at other times during the 
simulation. Chapter 4 will contain further discussion of these adjustments 
for the techniques considered during this research. 

3. 3 IFAM Logic Design 

The efficient implementation of the IFAM test bed on a computer 
with limited core storage such as the 200Kg UNIVAC 1108 was a significant 
problem. Both forward and adjoint data would be required during a single 
run. Additional core storage would be required by the importance function 
for the current mode and the flux estimator to generate the importance 
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function for the next mode. Both of these data arrays are functions of 
three parameters: position, energy, and direction. Additional storage 
area is required for manipulation of the raw importance data to prepare 
it for the next mode calculation. The core storage required for a typical 
problem with 40 importance regions, 20 energy groups and 18 angular bins 
is 14, 400 words for the importance function. The same number of words 
are needed for storage of the data which will be used to generate the 
opposite mode importance function. Additional storage is required for 
marginal probability distributions of the importance function. If scratch 
areas are set aside for data manipulation, then the core storage require- 
ments for IFAM - peculiar data exceeds 44, 000 words for either mode. 

The total for both forward and adjoint modes is greater than the 65, 536 
(200Kg) words of UNIVAC 1108 storage. The 44, 000 words does not include 
stor age for the instruction bank and for random walk, geometry, cross 
section, and analysis data. 

Obviously, techniques for reducing core storage requirements 
are required. Use of readily available techniques such as segmentation 
of che program provided some reduction, but not enough. Data could be 
stored on data files, but the increase in run time and input/output operations 
to retrieve and store data from data files during the random walk compu- 
tation was intolerable. It was obvious that all data needed for a given mode 
calculation must be in core. To meet this condition and to allow sufficient 
core storage for the additional data arrays, the logic flow of IFAM was 
designed as shown in Figure 5. This design is such that both segmentation 
and data storage or retrieval on data files can be utilized. In addition, 
only that portion of the input or generated data that is required during a 
given mode calculation resides in core. Scratch areas of core have been 
eliminated by dynamic allocation of data array storage for temporary data 
arrays. The result of this design is that the amount of core storage 
available for data is more than 18, 000 words greater than that for the 
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Figure 5. IFAM Logic Design 
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MORSE codes even though IFAM requires a larger instruction bank. Thus, 
iterative forward-adjoint Monte Carlo problems can be executed on the 
UNIVAC 1108 without removal of any of the MORSE code capabilities. 
Subsequent paragraphs will describe the IFAM logic design in more detail. 

As explained in Section 3. 1, the IFAM test bed is built around 
the MORSE code. The overall design of IFAM is based on the calculation 
of both forward and adjoint MORSE solutions alternately for a specified 
number of iterations. Before ea*. » mOue calculation, an importance function 
is computed and used during the mode calculation to alter the sampling 
distributions. This importance function is derived from the flux estimate 
stored during the Monte Carlo simulation of the opposite mode problem. 
Figure 5 is a simplified diagram of the logic flow of this design. Each 
process, input/output and decision operation have been nunbered for ease 
of reference in the following discussion. Some of these numbered operations 
are much more complex than others. For example, operation 14 consists 
on one line of FORTRAN code, while operation 2 and a subset of operation 
10 encompass a complete MORSE forward mode calculation. However, 
these operations do illustrate the major logic design of IFAM. 

The initial IFAM - peculiar operation occurs in the executive 
or controller subroutine for IFAM, MORSE A FA. Iteration control and 
source region information are obtained by input operation 1. The iteration 
control data includes the number of iterations, initial mode, the number 
of batches and the number of particles or adjunctons per batch for each 
iteration. Operation 2 is the input of the forward or F-mode data. This 
data includes the random walk, geometry, cross section and analysis data 
as described in Appendix B. During the F-mode data input operation, both 
the combinatorial geometry and the cross section data are temporily stored 
on data file 16. This data is read during the A -mode input operations. 

After all F-mode data has been read into the appropriate common areas. 
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the FAMS subroutine is called and the data in these common areas are 
written onto data file 17. This data file always contains the F-mode data 
which is transmitted back into the appropi ~*te core storage data array at 
the beginning of a F-mode calculation (represented by operation 14). 

Following the F-mode data input is the A -mode data input. The A -mode 
must follow the F-mode data input operation, but the initial random walk 
calculation is arbitrary. The A-mode data is read from cards (or card 
images) which are on the input file and data file 17. A-mode data is 
stored in the same areas of core storage as the F-mode data, and likewise, 
written onto the A - mode data file (16) by the FAMS subroutine. 

After all input data has been read and stored on the appropriate 
data file the iterative forward adjoint calculation begins. The initial 
mode (forward or adjoint) is determined during operation 8 and the data 
for that mode is retrieved from its data file. For example, if the initial 
mode is forward, then data file 17 is read and the F-mode data put into the 
core storage data arrays. This operation (15) is performed by the FAMS 
subroutine. FAMS also initiates operation 9, which consists of calling the 
FAIF subroutine for determining the importance functions and initializing 
certain variables and parameters. With these operations, IFAM is now 
ready to begin the random walk and analysis steps that constitute the irtajor 
part of the computation. The logic flow for this part of IFAM is very 
similar to the MORSE random walk, except for modifications to accommodate 
the distribution functions which have been altered by the importance function. 
In addition the data which will be used to generate the importance function 
for the opposite mode (i. e. , the track length or collision density) is also 
stored by angular bin. In MORSE, only energy group and importance region 
dependent quantities were stored. A fter each source selection and collision 
event, flux estimates are made to the point detectors using statistical 
estimation techniques. This data is cumulated and then output at the 
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completion of operation 10 for the initial mode of the first iteration, the 
equivalent of one complete (but short) MORSE run is terminated. 

Of course, IFAM does not terminate, but performs operation 11 
of Figure 5. This operation, executed in the FAMS subroutine, is 
the first step in generating the importance function for the next 
mode calculation. It is followed by an output of the initial mode 
data at the completion of the calculation to the appropriate data 
file. While many of the data arrays written on the data file are 
the same as the input data for that mode, some arrays will have 
been updated. These arrays include the geometry arrays, the 
scattering and track length counters, and the user (analysis) arrays. 
The data arrays at the completion of this initial mode are then stored 
on the initial mode data file, as indicated by operation 12. 

Since the above discussion was for the initial mode calculation 
for the first iteration, then decision operation 13 will be negative, 
and the mode will be changed to the next or opposite mode from the 
initial mode at operation 14. The data stored on the data file for this 
mode will be read and placed in core storage. For example, if the 
adjoint mode had been specified to be the initial mode, then at 
operation 15, the F-mode data is stored in c ire so that a forward 
Monte Carlo simulation can be performed by operation 10. Assuming 
that more than one iteration is requested, then the next time that 
operation 15 is executed, A-mode data will be put into core. In 
any case, data for the "non-initial’’ mode is now in core storage. 

The importance functions for this mode is calculated in the FAIF 
subroutine and is passed to the random walk routines by operation 9. 

At operation 10, another complete Monte Carlo calculation is 
performed, including the output of the results. This calculation 
is for the opposite mode from the Initial model, but uses data from 
the initial mode calculation in the importance function to alter 
the sampling distributions. 
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After operations 11 and 12, a check is made to determine 
if all iterations have been completed. If not, then the mode is 
changed again, this time from the "non-initial" to the initial mode. 

The initial mode data is stored in core, new importance functions 
are calculated from the previous run data, and the Monte Carlo 
simulation routines of operation 10 are executed. After operations 
11, 12 and 13 are performed, the "non-initial" mode calculations 
for the second iteration are executed. 

This procedure continues until both mode calculations are 
performed for all iterations. Then decision operation 13 is affirmative, 
and the end -of -problem results are printed. Note that the above 
logic design never requires both forward and adjoint data in core 
storage at the same time. In fact, the only significant amount of 
additional data is the array used to store the tlux estimator (WATK). 
Further reduction in core requirements is obtained by segmentation 
of the program into a main segment and four subsegments, with 
nearly 19,000 words saved. Appendix A contains a more detailed 
description of the segmentation of IFAM. Finally, the requirement 
for scratch area in core was overcome by storing data on the F- and 
A- mode data files until it could be placed in the blank common area 
used by the flux estimator without destroying useful data. This 
logic design resulted in a factor of 3 reduction of total core storage 
requirement, making possible the utilization of the IFAM test bed on 
the 65K word UNI VAC 1108. 

3. 4 Special IFAM Subroutines 

Twenty two of the subroutines contained in the MORSE code had 
to be modified in order to be compatible with the IFAM test bed. In 
some cases, the subroutine had to be completely revised for IFAM 
(e. g. NXTOOlA while other required only minor changes. In addition 
to these twenty two subroutines, three new subroutines were required. 
Two of these subroutines, FAMS and FAIF, are required to manipulate 
the data between the two mode calculations and to generate the import- 
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ance functions. Thethir'* cubroutine, ANGBIN, determines the angular 
bin for directional dependent storage or selects an outgoing angular 
direction isotropically over a specified angular bin. 

These subroutines have functions that are central to the design 
of IFAM, and are described in the following sections to provide a better 
understanding of the IFAM test bed. 

3.4,1 Subroutine FAMS 

Very early in the design of the IFAM test bed, it was recognized 
that a software routine would be required to handle the exchange of 
forward and adjoint data between mode calculations and control the 
calculation of the importance functions. Consideration was also given 
to using the same routine as the executive controller for IFAM, but 
further analysis revealed that this function could be handled easier by 
modifying the MORSE code controller (which is the MORSE subroutine). 
The modified MORSE subroutine, MORSE/I FA, is now used as the 
executive controller for IFAM, but it calls the FAMS subroutine for 
mode dat. exchange and importance function control. Almost all of 
the data used in FAMS is transmitted through the common areas. Only 
two parameters are passed through the subroutine agreement list, 
lADJ and MSR. IADJ specifies the current mode under consideration. 
MSr specifies whether the mode data is to he stored on or retrieves 
from the data flies. The value of MRS also determines what data 
manipulation will occur, as shown on Figure 6. 

During the read operation, the flux estimate from previous mode 
calculations is read into the FAI data array, which contains the import- 
ance function during the random walk. For the first iteration and 
initial mode, where no estimate exists, the FAI array is initialized 
to zero. The design of the FAMS subroutine is such that either the 
cumulative sum of all appropriate mode calculations or just the pre- 
vious mode results are retrieved. Next, the input data is read from 
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the mode data file and loaded into labelled and blank common areas of 
core storage. Only that data that is required for this particular mode- 
ls put into the core. This read operation performs the same function 
for 1FAM as the input operation does for the MORSE code. It is follow- 
ed by a call to the FAIF subroutine, where the importance functions 
are generated for the present mode calculation. The initialization of 
blank common storage areas and other variables complete the read 
option tasks of FAMS. The store option (MSR<0) tasks are performed 
at the conclusion of operation 10 on Figure 5, which include the random 
walk calculations for the current mode. The flux estimate determined 
during the random walk is either transferred to the FAI data array or 
it is summed to all previous estimates of the flux for this mode (for 
cumulative sum importance functions). In either case, this data is 
written onto the current mode data file along with the other labelled 
and blank common data arrays. The final flux estimate is retrieved 
from this data file during the next mode calculation, and used to generate 
the importance functions. 

FAMS will also print the contents of all the labelled common data 
arrays which it handles and of selected portions of blank common. This 
print is controlled by input data and provide a valuable aid in debugging 
operations for both coding and input data checkout,. 

3,4.2 Subroutine FAIF 

The generation of the importance functions which are used to 
alter the source, transport and collision distribution functions is 

performed in the FAIF subroutine. While the exact design of this sub- 
routine is dependent on the particular form of the importance function 
being tested, several general features of FAIF are common to most 
cases. These features are described in the context of the version of 
FAIF that normalizes the importance function array FAI and the input 
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energy importance array EPROB on a region basis. Justification of 
any particular method of generating the importance functions will not 
be discussed until Chapter 4. 

Figure 7 indicates the major features of the FAIF subroutine logic 
flow. During the initial mode of the first iteration, no flux estimate 
exists, so only the input energy importance array EPROB requires 
normalization. This data is then used to calculate the initial importance 
functions. The functions calculated are: 

1. RI - The region importance function 

2. FAI - The energy and angular importance function normalized 
for each region 

3. FAI + - The energy importance function for each region 

4. BFS - The biased source energy function. 

Details on how these functions are employed in IFAM are given 
in Chapter 4. However, it should be noted that FAi is generated by 
summing the FAI data array over all angular bins. The biased source 
energy function BFS is calculated by multiplying the input source energy 
function by FAI + for the source region and then normalizing. The above 
functions are calculated prior to each new mode calculation. 

Whenever the FAIF subroutine is called after the initial mode of 
the first iteration, the flux estimate from the opposite mode random 
walk must be normalized and stored in the correct angular bin for the 
current mode importance function. Storage in the correct angular bin 
is required because the flux estimate was stored in such a way that the 
particle or adjunction is traveling along its velocity vector instead of in 
the opposite direction as required by the importance function (explained 
in Section 2.4.3 and Appendix D). Correcting the angular direction or 
bin is handled by a transformation or permutation array defined in 
FAIF. The flux estimates are also divided by the region volume. It 
should be noted that the flux estimate has been stores in the FAI array 
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Figure 7. FAIF Subroutine Flow Diagram 
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by the FAMS subroutine. Also, all operations on the flux estimate, 
including the angular bin correction, are performed in a manner that 
requires no additional storage area. 

Once the flux estimate has been properly normalized and trans- 
formed, the final form of the importance function is calculated by a 
linear combination of the input energy importance function EPROB 
and the flux estimate. The final form of the region importance function 
HI is also calculated. RI is stored in the unused fission weight array 
FWLO of blank common. Next, the FAl^ and BFS arrays are calculated 
from the FAI array. Thus, the importance functions and, in one case, 
a biased distribution function (BFS), are computed by the FA1F sub- 
routine for the iterative forward-adjoint Monte Carlo random walk. 

3. 4. 3 Subroutine ANGB1N 

In order to satisfy the requirement for angular biasing with the 
importance function, the flux estimate taken during a given mode 
calculation mus^ be stored as a function of direction. The MORSE code 
provides for angular data at detector points, but not as part of the 
region and energy dependent random walk arrays. Therefore, it was 
necessary to set aside core storage for the region, energy and angular 
dependent flux estimate array. To handle the angular dependence, an 
angular bin structure was designed as shown in Figure 8. Also, the 
ANGBIN subroutine written to place each particle or adjuncton in the 
correct bin In addition, a software routine is required to select the 
specific direction of an emerging particle or adjuncton after an out- 
going angular bin has been selected. ANGBIN also performs this task. 
Figure 9 illustrates how these two tasks are accomplished. Discussion 
of the methodology used by ANGBIN and a description of the ANGBIN 
operations are given below. 

The angular bin structure was defined oy defining eighteen equal 
solid angle bins. These bins can be defined on a sphere (see Figure 8) 
as follows: 


3-26 



2 



Figure 8. A ngular Bin Structure 


$-27 





3-28 












Bin 1 - Polar cap about the *z - axis which is subtended by 

the cone whose polar angle has a cosine of 8/9. 

Bin 2 to 5 - Surface areas lying between the polar angles whose 
cosines are 8/9 and 4/9, each of which occupy one 
quadrant, with bin 2 in the first quadrant, bin 3 in 
the second, bin 4 in the third, and bin 5 in the fourth. 

Bin 6-9 - Surface areas lying between the polar angles whose 
cosines are 4/9 and 0 (the equator), each occupying 
one quadrant in sequential order. 

Bins 10-18-These bins lie in the -z hemisphere and were deter- 
mined by taking the absolute value of the z-axis dir- 
ection cosine, locating the -t-z bin, then determining 
the bin number by the formula: 

bin m - 19 - n, where n is the bin. 

This angular bin structure allows for the rapid determination of the 
bin from the direction cosines of the particle or adjuncton path. Each 
angular bin has a solid angle of one -eighteenth of the total solid angle 

(i.e. 2 ff ). Once the polar region of the direction has been determi- 

T' 

ned, the quadrant is found by checking the sign of the x- and y- 
direction cosine. This eliminates the necessity for any complex arith- 
metic calculations. Figure 9 shows the logic flow for determining the 
angular bin from the direction cosines. 

ANGBIN also selects the direction cosine uniformily over the 
input angular bin whenever the option variable, IOP. is input as zero. 
Since each angular bin is limited in its polar and azimuthal angle extent, 
then the cosine of the polar angle and tb P azimuthal angle are selected 
uniformily over the limits of the input bin. This can be done quite simply 
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by adding ihe product of an uniformity distributed random number times 
the angular interval to the lower limit. If the lower azimuthal angle of 
the input bin is ft ^ , and the bin width is & ft^ , then the azimuthal angle 
ft is calculated by: 


ft = ft + R . £ o , (3 - 6 ) 

l n l 

where R^ is the random number. Since the polar angle is cosine dis- 
tributed, then the cosine of the polar angle is uniformity distributed. 
Thus, the z - axis direction cosine C can be computed from knowing 
the cosine of the lower limit C ; and the cosine width AC^. Selecting 
another uniformily distributed random number R n , C is calculated as 
shown below: 

C = C < R . A C . (3 - 7 ) 

l n 

The direction cosines for the x - and y - axis can now be calculated 
from the trigometric relations: 

A = cos ft , (3 - 8 ) 

B -- sin ft \ll - c 2 , (3 - 9 ) 

/ 2 

where Vl - c is the sin of the polar angle. In this manner, 
ANGBIN determines the direction cosine (A, B, C) of a direction 
uniformily distributed over the input angular bin. 
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4. 


METHODS DEVELOPMENT 


The theoretical discussion of the iterative forward-adjoint 
Monte Carlo method in Chapter 2 demonstrated the validity of the 
method for biasing the sampling distributions in continuous phase 
space. However, the solution of most practical radiation transport 
problems requires the use of high speed computers and the reduction 
of continuous phase space to finite phase cells in energy, position and 
direction. Chapter 3 has already introduced the phase cell structure 
which is employed in IFAM, but it did not answer the important question 
of the validity of the finite c' Hs approach. This question can only be 
answered by the development of methods for applying the information 
obtained during the opposite mode calculation and then testing these 
methods on real radiation transport problems. Methods development 
will be addressed in this chapter, and the result of applying these 
methods to a radiation transport problem is given in Chapter ?. 

The central problem in the development of methods for biasing 
of our sampling distributions is to defi the importance function. 

This problem is discussed in Section 4.1, where the form, format, 
and restrictions which were placed on the importance functions are 
described in detail. Application of the importance function to the 
sampling distributions is discussed in the order that IFAM employs 
them. Section 4. 2 describes the biasing on tne source distribution 
in both energy and direction. The transport of the particle or adjuncts 
is next in IFAM, so Section 4. 3 defines various techniques considered 
for the bia ing of the transport kernel. The mathematical justification 
of these techniques is given and inherent problems and restrictions 
are discussed. The major problems encountered during this research 
were caused by the use of importance sampling techniques on the 
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transport kernel. However, since many important radiation transport 
problems (e. g. , any involving deep penetrations) require altered 
transport kernels, considerable effort was spent in resolving these 
difficulties. 

Following selection of the coll ion ite with the altered trans- 
port kernel, the outgoing energy ana direction must be selected from 
the collision kernel. Methods were also developed for altering the 
sampling from the collision kernel by using the importance function. 
These methods are discussed in Section 4.4. The first four sections 
emphasize the development of importance sampling methods. Section 
4. 5, which is independent of importance sampling, contains a brief 
description of a technique for avoiding the infinite variance problem 
when estimating the response at a [xnnt detector with the statistical 
estimation or expected value metiiod. This technique was required 
because the detector is contained in scattering region, thus allowing 
the possibility of an "infinite’’ estimate. 

4. 1 The Importance Function 

The theoretical application of an importance function for altering 
or biasing the sampling distributions in a Monte Carlo radiation 
transport calculation were discussed in Section 2. 3. In Section 3.4, 
the computer routines which generate the importance function were 
described in sufficient detail to explain the logical flow of the 1FAM 
test bed. In this section, the analytical and computational rationale 
for the specific forms of the importance function will be given. As 
stated in Section 2. 3. the importance function was used to alter the 
sampling distributions, but the bias introduced was removed by correct- 
ing the weight (e, g. , the source particle is selected b om the altered 
source distribution S'(P), but the weight assigned to that particle is 
the ratio of the natural or true distribution, S(P), to S’(P)). 
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The initial consideration in developing the importance function 
was what quantity should be stored during a given mode (forward or 
adjoint) calculation lor defining the importance function for the opposite 
mode. The theoretical considerations o f Chapter 2 indicated that a 
proper quantity would be the adjoint to the integral equation being 
simulated by the Monte Carlo procedure. Analysis revealed the f-ict 
that the adjoint angular flux, y * (x , u: ) or <i>* (x , ~), is adjoint .o the 
emergent particle density equation (2-92). The forward angular flux 
c p (x , w), is the quantity "adjoint’’ to the integral equation (2-117) 
similated during the adjoint mode ( 1 . e. . the "adjoint equation is a 
"forward’' quantity, the angular flux). Table 1 of Chapter 3 shows 
that each of these two quantities can be estimated by two different 
techniques. The first technique is based on the collision density 
t.. e. , o = £ t p). The forward or adjoint flux is estimated by dividing 
the statistical weight at a collision point by the group cross section. 

The second technique is based on the track length in a region. Since 
the average flux in a given region is equal to the total track length of 
all particles crossing the region divided by the region volume, then 
the flux (forward or adjoint) can b< estimated during the Monte Carlo 
simulation by dividing the product of the statistical weight times the 
path length in a given region by the region volume. Both of these 
techniques have the attribute of requiring the same code design for 
both the forward and adjoint mode calculations. 

For a given test run, only one of these techniques can be used 
to estimate the forward or adjoint angular flux. Therefore, it is not 
necessary to store both the collision density and track length flux 
estimators, and core storage is set aside for only one of these f sti mates. 
The track length estimator has been examined most thoroughly, since 
it provides an estimate of the flux in a region whenever the particle's 
path intercepts the region, although no collision occurs there. The 
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:rack length estimates are stored in the I FA version of the NXTCOI 
subroutine. The statistical weight used to multiply the tt.uk length 
in each region is computed at the next collision point or at the escape 
point. Track length data is stored in NXTCOL/IFA until this weight 
is determined, and then a running sum of the product of the weight and 
track length is stored as a function of energy, region and direction. 
Division by the region volume is performed for the sum of these estimates 
in the FAIF subroutine. 

Storage of the collision density estimator occurs after the trans- 
port step (CALL 1IATCOL) and before the collision step (CALL COLISN) 
n the MORSE subroutine. Like the track length estimator, a running 
sum of the statistical weight is stored for the energy group, region 
and angular bin of the particle at the collision point. This sum is 
divided by the group cross section to produce the proper angular flux 
estimate. Because of the similarities between these calculations and 
those for the MORSE scattering counters, this data has been put in 
blank common immediatel / following the energy and region dependent 
track length array. 

At the beginning of the subsequent mode calculation, this data is 
summed with the track length or collision density estimates from 
previous IFAM interations. Provision has been made to use only the 
previous mode data for the estimates, but the statistics are very poor. 
This practice was soon discounted as a reasonable method of determining 
the track length or collision density estimates for the importance 
function. 

Although the theoretical analysis indicated that the adjoint to 
the integral equation being simulated was a proper choice for the 
importance function, computational considerations place certain restric- 
tions on this choice. The initial restriction is due to the fact that it 
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is highly unlikely for each finite cell in phase space to have an accurate 
estimate of the forward or adjoint flux. For example, consider the 
test problem, which has 15 energy groups, 42 importance regions, and 
18 angular bins, or a total of 11, 340 phase cells. Even if the histories 
were uniformily distributed over these cells (of course, they aren't), 
ten percent accuracy would require over a million estimates. This 
is an impractical requirement, so a procedure developed to use 
already available importance information in addition to the Monte Carlo 
"adjoint function" estimate. This procedure can be represented by the 
following equation: 


I(i, j,k) = a n J(i, j, k) + (1 - a n ) C(i, j, k). 


(4-1) 


where 


(i, j,k) 
I(i, j, k) 


J(i, j, k) 
C(i,j,k) 


indices of the phase cell for the enerrv group, region, 
and angular bin, respectively, 
importance function for the (i, j, k) phase cell, 
a parameter which is dependent on the number of 
iterations or histories that have been executed 
(0 < a n • 1), 

normalized angular flux estimate from the opposite 
mode calculations for the (i,j,k) phase cell, 
normalized input value for the (i, j,k) phase cell. 


Equation 4-1 allows the IFAM test bed to increase the reliance 

on the "ad’cir.t" flux (whicn is the forward anrular flux for altering 

the distributions for the adjoint integral equation simulation) as the 

statistical accuracy of that estimate increase. This is done by 

increasing a . One method that has b en used is to define a by 
n n 

n/(n + 1), where n is the number of iterations. Another possibility 
is to pick some minimum number of histories, say m. then set to 
n/(n + m), where n, is the total number of histories completed for the 
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specified mode. The importance function thus approaches the .If i , j,k) 
value as the number of histories increase. Also, r equiring that C(i,j,k) 
must be positive ( 0) and that must be less than 1 assumes that 

no cell in phase space will have zero importance. This is a potential 
danger when enly .J(i, j,k) is used in the importance function. One very 
good choice of C(i,j.k) is the adjoint flux from discrete ordinates 
calculation where the finite cells in phase sitace used in the discrete 
ordinates and IFAM are matched as closely as possible. 

The importance function defined by Equation 4-1 provides the 
basis for altering the Monte Carlo sampling distributions in IFAM. 

As discussed above, there exist several options for determining each 
of the major parameters in the equation (a^, .1(1, j, k) and C(i, i, k ) ). 

The result of chosing a representative range is given in the next chapter. 
In addition, the most useful form of the importance function depends 
on the particular sampling distribution being altered. The next three 
sections discuss the forms in which I(i, j,k) is appL d. 

4. 2 Source Biasing 

Selection of a source particle usually requires that an energy 
group, spatial position, direction, and ,.ge bo selected from the distri- 
bution functions. For :he IFAM test bed, th-- age is assumed to be 
zero and the spatial position assigned at the input source point. Then 
an energy group a d direction is selected from the altered distribution 
functions. While selection of an energy group and angular bin is 
possible from a joint probability distribution functie , the method used 
in IFAM is consistent with the MORSE procedure cf selecting the energy 
group from the marginal distribution function, then selecting the angular 
bin from the conditional probability function given in the above energy 
group. The specific direction is defined by selecting the direction 
cosines so that thi direction vector is uniform over the previe sly 



chosen angular bin. Details of this methodology are discussed below. 

Consider first the selection of an energy group. The natural 
energy dependent source distribution is one of the default input options 
for IFAM, To generate the altered energy group distribution function, 
the energy and region dependent marginal distribution function is 
computed as shown below: 

r(i,j,k) = I(i,j,k)/V V i( i)jtk ) (4-2) 

i k 

where I(i, j,k) is defined by Equation 4-1. This new importance function, 
I'(i, j.k) is now normalized on a region basis. Next, the marg nal 
distribution for the energy groups and importance regions is computed 
from I* (i > jj k): 

I' k (i,j) (4-3) 

k 

The final step in defining the altered energy dependent source 
distribution, S (i), is: 

S(i)= S(i) • I R (i,j s )A[S(i) • I k (i,l s )], (4-4) 

K 

where S(i) is the input natural distribution and i is the importance 

D 

region which contains the source. Note that the distribution has 
been normalized, so that the initial statistical weight for a source 
energy group of i g is given by: 

W g = S(i S (i g ) . ( 4 - 4 ) 

Only the source direction remains to be determined. As in the 
source energy selection, the importance function for the region contain- 
ing the source point is used to supply angular biasing information. 

Since the energy group and region are known, then the proper importance 

function to be used is I' (i , j , k). Defining D(i , j , k) as the angular 

s s s s 

bin distribution function, then the angular bin is chosen from: 
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(4-5) 


D (i g , j s , kl D ( i g , j g , k) • 1' (i s , j g , k) 

- , L. k)I’ (i j k)] . 

L ^ ^ » a j 

After the angular bin, k , has been randomly chosen fr >m (4-5), the 

s 

direction cosines from the direction vector are determined in the 
ANGBIN subroutine, as discussed in Section 3.4. However, the use 
Oi D (i g , j g , k) instead of D (i g , j g , k) for the selection of the outgoing 
angulai' bin means that the statistical weight must be corrected. The 
final statistical weight for the source is given by: 


W 1 • D(i ,i , k I ^ (i , j , k ) . 

s s s’ s s J s s 


For the special case where the source is isotropic. Equation 
(4-5) and (4-6) can be written in a much simplier form: 


D<i S ’V k) = ■’ <vv k)/ V> S ') S > ■ 

w f = vv l /fk • i; u , j )1 , 

c c I m o v Ir 3 


(4-7) 

(4-8) 


where I’ <i , j ) is defined by Equation (4-3) and k is the total 

K S S IY13.X 

number of angular bins. The two equations above were used for the 


test problem. 


4. 3 The Altered Transport Kernel 

The r„ ost difficult of the regular Monte Carlo sampling distri- 
butions to alter with the ge *rated importance function is the transport 
kernel. The reason is that the transport kernel is a continuous 
exponential function, whereas the importance function is a set of 
numerical values, one for each finite cell in phas? space. The 
importance function is easily applied to the discrete sample space of 
the source and collision kernel, but this is not the case for the 
transport kernel. To better understand the mathematical foundation, 
consider the following form of Equati n (2-62): 
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X'(x , E) = S (x, E) + / C (E | E'; x) T (x |x’; E’) Y(x’, E')dx'dE\ 

(4-9) 

All quantities are defined the same as m Chapter 2. 

Now, as in Section 2. 3, multiply the equation abo’ e by the 
importance function, represented by I (x, E). Letting: 

V(x, E) I (x, E ) y(x . E) , (4-10) 

and S (x, E) I (x, E) S (x, E), (4-11) 

then: 

y(x, E) = S (x, E) + fc (E | E' ; x) T (x|?:E'' 

[i (x, i)A (x’, E’)J V (x ' , E’) dx f dE’. (4-12) 

Consider the following part of the integral term of (4-12) 

C(E | E' ; x ) T (x | x'; E') [i (x. E)/I (x',E’)J 
= [i (x, E)/I (x, E’)] C (E | E'; x) 

• [i (x, E’)/ * (x\ E’)J T (x I x E’ ) . (4-13) 

Inspection of the right hand side of the above equation shows that the 
importance function ratio has been separated into the product of two 
ratios, the first of which seems to apply naturally to altering the 
collision kernel and the second to the transport kernel. For example, 
the ratio multiplying the transport kernel is just the ratio of the 
importance functions at the beginning of fhe transport step (x’) to that 
at the terminal point (x). The directed energy (E’) is that which 
describes the energy and direction during the transport step, as 
would be expected from physical considerations. 

To use the altered transport kernel, some method for sampling 
must be devised., This is a non-trivial problem, since the representation 


4-9 



of the kernel in thk Vlonte Carlo simulation for the discrete phase 
is: 

n-1 

T ( 0- R)dR - a n .Ij) - n ’ xp -[N_^ j (R ! - R j _ J ) 

* V (R - R dR 14-14) 

n n-1 J 

where T is the importance function value in tae i -th region which the 
path of tne p. rticle intersects, with the collision assumed to occur 
in the n-th region. The total group cross section for region i is 
represented bv u., and R. is tnc distance to the i-th region boundary. 
The distance from the source point to the collision point is R, and R q 
is defined as 0. If the collision occurs in the same region as the 
starting point, then the summation term (from j - 1 to 0) is just 0. 
Note that ii the importance function ratio in the right hand side of 
Equation (4-14) is dropped, then the remaining statement is jvst the 
natural transixirt kernel, T (0 - R)dR. 

The problem with sampling from Equation (4-14) is that the 
kernel is not in general normalized, that is: 

oo 

f T (0 - R)dP * 1 (4-15) 

*0 

This can be seen from the following representation: 

i-i 

V u. (R. - R. „ 

*— J ' J-l) 

H 

(4-16) 

J T ( 0 - R) dR - i' 1 Jlj - (I 2 - Ij) exp -(p^) 

2 

+ (I 3 - I 2 ) exp -(^ u. (R. - R ^ _ J ) ) e J (4-17) 


R. 


f T 1 1 

“'O 


(n-R)riR =V f (!./!,) fi. exp 

m-L ' 1 1 


+ u. ( R - R. .) 

l l-l 


rii 

dR 


i-i 


4-10 



i 

M j (R r R i-i ) J * 

(4-18) 

Note that for the special case where all L’s are equal, then the 
altered kernel reduces to the natural kernel and is, of course, 
normalized* 

By forcing each path to an escape boundary and determining the 
value of Equation (4-18), a normalized kernel can be constructed, but 
the additional computation time will be excessively large for problems 
with many regions. If this procedure is used, then the new kernel is 
defined by: 

* tt 

T(0-R)dR = T(0-R)dR/£ T (O-R)dR, (4-19) 

where the integral must be evaluated for every transport step. 

The usual procedure employed in the alteration or biasing of the 
transport kernel is to multiply the total group cross section term, 

Mp by a fixed value. This technique can be extended to multiplication 
by a variable which depends on the position. For the problems of 
interest, the assumption of a region dependent parameter is valid, 
as shown below. Let: 

(4-20) 

where is the region dependent parameter for a given path. This 
means that Cj can be energy, direction and even initial position 
dependent. Defining the altered transport kernel by: 

T ( 0 - R) dR * u j exp -[£ |*|(Rj - R^) + (R - R ul )] dR, (4-21) 

produces a normalization factor of: 


or J T(0 - R) dR = 1 + I" 1 . £(I - I ) exp -[£ 

0 i=l L j=l 


^ 


4-11 



OB 1 — 1 

f T(0.R)dR=X/ C.4j exp -[£ Cf.Ol.-R._J 
° i=1 Ri-1 j=l 

+ CjU. (R-R. X )J dR 


(4-22) 


-Z !»•>-[£ c i“i (R r R j-i ) ] I b i“i exp -[ c i“ 1 ( R i- R 1 -i>] dR l’( 4 - 23) 
i:} R i-i 

=X e *p -fz c i‘‘) <R i' R i-i ) ] • 1 • exp " [ c i , 1 i <R r R i-i ) ]i ’ <4_M) 

i-1 |=l J 1 L 

=z | e *i> -[z c i M j (R i - R j-i ) ] - exp -[z c j ,, i (R j- R j-i>] | • <4 - 25 ' 

Expansion of the above equation reveals that the negative exponential 
term for i = n is the same as the positive exponential term for i = n + 1. 
Thus, each of these pairs cancel, leaving oi ly he positive term for 
i = 1. Since the summation from j = 1 to O’, was previous shown to be 0, 
then the normalization factor is the exponential of 0 or: 

«l 

^T(0-R)dR = l. (4-26) 

This proves that any altered transport kernel which is based on 
replacing the total group cross section with the product of a region 
dependent parameter C and is normalized. An example of (Z is 
the reciprocal of the constant path length stretching parameter BIAS 
used in MORSE code. (Note that although BIAS is computed as a 
region, energy and direction dependent parameter, only the one value 
for the starting point of the path is used to alter the cross section 
for all subsequent regions intersected by that path. ) Another example 
of C, is the region importance function Ij for a specified path (and 
hence energy group and angular bin). Usually ^ is put through some 
normalization process, such as: 

C, ( 4 - 27 ) 
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(4-28) 


°r C t = 1;^ • 

j J 

The summation of j is for the regions along the projected path of the 
particle or adjuncton, and x j could be the mean free path through 
region j. 

Two methods have been examined for the implementation of the 
importance function in the altered transport kernel. Both methods 
use the concept of total cross section alteration by Cj, since the 
normalization factor for the altered kernel is unity. The first method 
will be described without including the path length stretching capabilities 
available in the MORSE code. However, these capabilities are included 
in the methods implemented in I FA M and their addition to the method 
will be discussed later. The notation will continue to denote region 
dependence only, since the energy and direction are fixed during the 
transport step. The major steps of the first method are outlined and 
described below: 

1. Select the number of mean free paths, MFP, from an 
exponential distribution. Both the natural and altered transport kernels 
are exponentially distributed, as indicated by Equation (4-21), where 

is unity in every region for the natural distribution. For this method, 
our next collision point will occur at MFP mean free paths along the 
path as determined by the altered cross sections. The collision point 
will, in general, be different from the point at which MFP mean free 
path is readied by using the real cross section along the path. 

2. Determine the search length, ETA, from a preselected input 
value. The purpose of ETA is to define the number of mean free 
paths along the projected path at which the importance function Ij, 
cross section and path length ^ data will be found. Beyond this 
search length, the cross sections will not be altered. ETA is 
arbitrary, but should be sufficiently large to assure that more 
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than one region is intersected before ETA mean free paths are 
reached. ETA should not be so large that the probability of reaching 
ETA mean free paths is very small (e. g. , 10 mean free paths). 

If MFP is greater than ETA, then a normal transport game will 
be played. 

3. Using the combinatorial geometry routines of IFAM, step 
through the geometry, determining for each region which the path 
intersects (out to ETA mean free paths): 

Ij - . importance function value for the i-th region intersected 
by the path (i = 1 is the starting point region) 
t^: track length through the i-th region 

U { ‘. group cross section for the i-th region 
\y mean free path through the i-th region (Xj=*ijtj) . 

The 1^, and Xj data are stored for future use. The region cross 
section is not stored since and t^ define u^, and is not needed 
in the computation technique . A: n mathematical convenience, 
the region containing the end of the search length ETA will be assigned 
two index numbers. That portion of thr path that is betwe en the 
starting point and ETA will be denoted l / m, where m is the total 
nr* tuber of regions intersected up to ETA mean free paths. That 
portion of the region beyond ETA will be m + 1, and that part of the 
path will not require an altei ed cross section. 

4. Calculate the altered mean free path through each region 
based on the constraint that the sum of the altered mean free paths, 

XJ, equals the sum of the real mean free paths, Xj. That is: 

z x i - l>i ’ <*-»> 

i=l i=l 
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(4-30) 


where xj is defined as: 

X 1 C X i x i ’ 

The constant of proportionality, C, for Equation (4-30) can be 
calculated from Equation (4-29): 

• (4 " S1) 

i=l i=l 

Inspection of the above equations reveals the feet that the altered 
cross section that corresponds to Xj can be written in the form of 
Equation (4-28), where: 

Cj = CIj . (4-32) 

Therefore, the altered transport kernel has the same form as Equation 
(4-21) and hence is properly normalized to unity. 


5. Calculate the distance along the path at which MFP mean 
free paths have been traversed as determined by the altered mean 
free paths in each region. This step is equivalent to selecting the 
distance R from the altered transport kernel. The computational 
procedure is shown below: 


a) Define the i-th region by: 

1-1 1 

£x[ s MFP < X \[ . (4-33) 

b) Determine the distance R from the starting point by: 

l-l l-l 

R “Z l i ♦ l i • < MFP - £ Xp/ X’ . (4-34) 

i*l t»l 

The calculation of R consists of accumulating the distance up to the i-th 


region and then adding that fraction of the distance that R goes into 
he i-th region. The collision point, r, in three dimensional space is 
calculated quite easily since die starting point, P , and die direction 
cosines of the path are known. 
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6. Since the above procedure produced a collision point different 
from that which would be obtained in selecting from the natural trans- 
port kernel, the statistical weight will be corrected to remove the 


bias. The weight correction is just the ratio of the natural to the 

altered transport kernel evaluated at the selected distance (R). 

1-1 

w _ T (0 - R) = 6XP ~[S X i * (R ~ ^-l*] , (4-34) 

T(0 - R) «' 4 exp -['f XJ + M' (R - R M) ] 

J=1 J 

where: . i 

"*-» * &•» • 

Since the ratio of u to is the same as X to x' , and the exponential 
factor in the denominator is MFP, then the weight correction is just: 



x i 

rr- exp (MFP - CMFP) . 
a j l 


(4-35) 


The term CMFP is just the true mean free path from the starting point 
to the collision point. It is computed in the same manner as R in 
Equation (4-34), with the t i and t £ replaced by X. and x lf respectively. 

The above procedure should produce a unbiased game, since the 
bias introduced by altering the transport kernel (Steps 1-4) are 
compensated by correcting the statistical weight. Additional changes 
have been made to the procedure to allow for the regular path length 
stretching. The technique implemented requires only that the search 
length ETA be multiplied by the bias term (BIAS) in step 2 and then 
Equation (4-31) which defines the constant of proportionality most 
include a division by BIAS. NO other changes are necessary to 
implement the use of MORSE path length stretching. Some 
logic was added to handle the special cases where the seerch leegth 
is contained in one region and where an escape occurs. The co mputatio nal 
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requirement for a single region path are much less than lor the above 
procedure. The logic for the escape condition is such that the above 
procedure is used unless the escape would have occurred in a regular 
MORSE run. 


The second method which was investigated for application of the 
importance function to the attempts to determine a pseudo-cross 
section which produces the same transport kernel as given by 
Equation (4-14). This was done by equating (4-14) to (4-21) and 
then solving for C^. That is, set: 


CjUj exp -[icj-VYV'] - IjAj) Uj exp-^u^-Rj j)] . (4-38) 


This equation can be written in the following form: 


C i = W R exp * [ (1 " C j ) u j (R i ' R j.i>] • 


Equation (4-37) is transcendental, but solution by trial and error 
seemed feasible. It can be shown that: 


(4-37) 


C i = VW C i-1 erp ‘ [ (1 " C i ) p i (R i “ “i-i^ «-») 

and C x = 1. (4-») 

Thus, it Is possible to solve (4-38) recursively. However, 
not all combinations of (IjA^) and u^Rj - R^) allow a real 
solution, so the use of this procedure must be sever ly restricted. 


A third procedure which may be impl e mented In future effort is 
b a sed on Equation (4-14) directly. As discussed earlier, this will 
require the generation of the normalisation factor. This can be done 
by forcing the path to boundary of the geometry model, collecting the 
track length and importance data, then evaluating Equation (4-17). 

It can also be done by projecting the path through a specified number 


4-17 



of regions and setting the importance ratio (I. Aj) to 1 for subsequent 
regions along the p&th. To preserve trie same escape probability, 
the importance ra!io for each of the specified number of regions 
may be renormal, zed. 

4.4 The Altered Collision Kernel 

The collision step takes the incoming energy and direction 
parameters at the collision point selected by the transport step and 
generates the outgoing energy and direction parameters. The major 
effort in altering the collision kernel with the importance function 
during this research has been to bias the selection of the energy 
group toward the more important groups. Equation (4*13) indicates 
that a proper importance function is just the ratio of the importance, 
I(X, E ), for the possible outgoing energies and directioi i to the 
incoming importance, I (X, E’). Representing the kernel in the finite 
phase cell induces format yields: 

[l (i c , j, k)A(i c , j c , k c )] . C(i c ; j e - j, k c - k) (4-40) 

as the altered kernel, where the subscript c on the indices denotes 
the incoming collision parameters. 

To select from the altered kernel, it must be normalised at each 
scattering event. This is done by adjusting the statistical weight to 
account for the absorption probability as discussed in Section 3. 2. 
Normalization to account for the importance function in the altered 
collision kernel is performed as shown below: 

C E (l c , J c - j) = N' 1 . Ijj (i c , i) . P <l c> ) c - j) (4-41) 

where I! (i , j) is defined by Equation (4-3). P(i , j - j) is the 

K C v v 

natural probability of scattering from group j into group j. 

v 
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(4-42) 


Ng is the normalization factor, given by: 

NDS 

N E = He ^c’ ^ P ^ ’ 

where NEK is the number of down scattering groups. This method 
was easily adapted to the available MORSE coding with the energy 
importance data being replaced by 1^ (i, j). 

Methods to alter the selection of the direction have not been 
implemented for the collision step. Work by other researchers 
have revealed that the technique of sampling from che importance 
fun ion suffers from the generation of negative weights (Reference 27). 
This is due to the fact that the reconstruction of the natural scattering 
distribution from the Legrendre moments results in negative values. 

For some problems, as many as 20 percent of the weights may be 
negative after correcting for sampling from the importance function. 
These researchers are examining the feasibility of altering Just 
the probability of scattering into the discrete angles used in MORSE. 

This technique avoids the negative weight problem, and preliminary 
calculations indicate a variance reduction for a highly directional 
problem. 

4. 5 Exclusion Volume Contribution 

The methods developed in this research are to be applicable to 

the difficult point source - point detector problem. This class of 

transport problems usually uses "statistical” or 'test flight" estimations 

to evaluate the detector response. However, tl tese problems suffer 

from die so called "infinite variance" problems whenever ti*c de te ctor 

is located in a scattering region. The infinite variance problem is 

.2 

due to the fact that the response estimator has an R term, where R 

is the collision to detector point separation. For collision very near 

-2 

the detector point, R is very large, and hence the cootrlbotion from 
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those collisions dominate the response estimate, and produce a much 
larger response than the actual value. One method of overcoming 
this problem is to place an exclusion region about the detector point. 
The response contribution from any collision occurring in the region 
is neglected. This procedure obviously produces low response 
estimates. 

To prevent the "infinite variance" problem and at the same time, 
include the contributions from collisions near the detector, the 
following technique was developed. Assume that che collision density 
within the exclusion region is uniform, then a collision of point r 
within the region could have occurred with equal probability at any 
point in the region. Thus, a collision with incoming energy and 
direction which has a statistical weight of will have the 
following flux contribution: 

W i 

<*(E)> = ^±-Jfu s (E i )/n(E i )HE i ~E,w i ~Zi) 


• exp (r’) dr'j /r 2 dw dr , (4-43) 

where: V = volume of the exclusion region 

CAC 

and P g (E { ) - scattering cross section for energy E^. 

If we make the exclusion region a sphe re of radius R, then: 


<*(E)> = W l 


exc 


. |~ P(E i - E)j J \ 2 exp -[jj*(R‘)dR* ]/r 2 


* /p(w j • w)dw dR . (4-44) 

The term in brackets outside of the integral is j jst the poet collision 
weight W , and the integral over the direction ta is unity. The 

B 

integral over R can be evaluated if a p (E) term is introduced. 
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This yields: 

<*(E» = [w g /(u ( E ) . V ex j\ . [l - exp (m(EVR)] . (4-45) 

This contribution is computed for each collision in the exclusion 
region and added to the total flux calculation so that the response 
due to these collisions can be estimated. 



5. 


RESULTS AND CONCLUSIONS 


The purpose of this research effort was to develop a computer test 
bed for examining the iterative forward-adjoint Monte Carlo method and 
to evaluate that method. Chapters 2, 3 and 4 contained a comprehensive 
discussion of the theoretical basis for the forward-adjoint Monte Carlo 
method, the development of the IFAM test bed and the variance reduction 
techniques that have been implemented in IFAM. In this chapter, some 
of the initial results will be discussed. Due to the large amount of effort 
required to develop, code, and debug the IFAM test bed, valid results were 
not obtained until the last month of this effort. Thus, the analysis of the 
iterative forward-adjoint method has been restricted to a few simple cases 
in order to meet contract delivery requirements. It is the intent of the 
author to continue the analysis of this method in much greater detail and 
for a large variety of problems. 

The results discussed here are for IFAM runs of three iterations 
with various importance biasing methods in effect. These methods include 
source direction biasing, source energy biasing and transport kernel biasing. 
For comparison, a run with no biasing and another run with all of the biasing 
methods discussed in Chapter 4 were made. The problem considered is a 
cylinder of air with a point fission source and point detector. The detector 
response function is the Henderson Tissue dose. Details on this problem 
are given in Appendix C. One difference that exists in the output in Appendix 
C and the problem results given here are the dimensions of the air cylinder. 
These results are for a source -detector separation of 300 meters. The reason 
for selecting this particular problem is the existence of extensive analytical 
and experimental results that can be used to validate the IFAM results (see 
References 28-30). 

Table II contains a summary of the IFAM results for the five runs. 

The uncollided dose is denoted by D . D is the contribution of collisions 

u’ exc 
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in the exclusion volume (see Section 4. 5) to the total dose estimate, and 
the total dose is denoted by D^,. The symbol FSD following each dose 
value stands for the fractional standard deviation, which is defined as 
the square root of the sample variance of the mean divided by mean value 
of the dose estimate. This quantity is the statistical estimate calculated 
by the MORSE and IFAM codes. The UNIVAC 1108 run time is also 
given, as well as the number of scattering events for each mode. Two 
estimates of the "goodness" of the variance reduction method are also 
presented. The efficiency is defined as the reciprocals of the square of 
the fractional standard deviation times the relative run time (which is 0. 1 
times the run time). The use of the relative run time was for convenience 
in plotting the results, and since the efficiency is a relative measure of 
goodness instead of an absolute measure, it has no effect on the conclusions. 
The second measure is denoted as the "% Error" and is based on the 
Henderson tissue dose calculated by the two-dimensional DOT -HI Discrete 
Ordinates computer code for the cylinder of air. That is: 

% Error = 100 • (D g - D T )/D g , (5-1) 

where D is the DOT -HI calculated value of 3. 408 X lCf Rads. Note that 

B 

this measure does not include any consideration for the computational effort 
required (such as the run time). 

The arrangement of 'Table II is such that both mode calculation results 

for a given iteration are adjacent and are in the same row as the results for 

the other cases for that iteration. The headings A and F denote the adjoint 

and forward modes, respectively. The total dose estimate, D^,, the total 

dose estimate from collisions outside the exclusion volume, D_ - D , and 

a me 

the efficiency has been plotted as a function of the lteration/mode sequence 
for the five test cases in Figure 10. 

Inspection of Tuble n and Figure 10 indicate that die utility of the 
iterative forward-adjoint method for the air cylinder problem is questionable. 
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TABLE H 


AIR CYLINDER CALCULATIONS 






Source Direction 

Source Energy 

Transport Kernel 





Mo Binning 

Binning 

Binning 

Binning 

AUBI 

Inning 



A 

F 

A 

r 

A 

F 

A 

F 

A 

r 



0.211 

0.214 

0.211 

0. 214 

0.211 

0. 214 

0. 210 

0. 213 

0.311 

0.309 


ISO 

.013 

.011 

.013 

.022 

.013 

.014 

.014 

.012 

.013 

.030 


D 

0.993 

1.391 

0.993 

1.989 

0.953 

1.616 

1.219 

1.197 

0.993 

3.130 


•“ no 

.921 

.156 

.927 

.213 

.827 

.223 

.930 

.397 

.on 

.744 

1 

V*** 

3. 399 

3.013 

2.399 

3.094 

2.359 

2.043 

2.939 

3.137 

2.399 

3.439 

I 

n r 

3.313 

3.393 

3. 312 

4.073 

3.312 

3.659 

3.997 

3.338 

3.313 

1.579 

2 

FSD 

.119 

.014 

.199 

.120 

.189 

.116 

.168 

.133 

.199 

.493 

3 

Ran Tim* 

400 

1030 

330 

970 

390 

1020 

390 

1030 

390 

900 


Scattering 

SUIT 

77131 

21937 

84903 

21937 

77135 

21994 

79910 

31937 

74930 


ItBclener 

0.999 

1.773 

0.949 

0.719 

0.717 

0.729 

J. 906 

0.843 

0.711 

9.933 


% Irror 

3.919 

1.330 

3.919 

-19.424 

3.918 

-7.335 

-12.499 

3.113 

3.919 

•83.973 



0.309 

0.313 

0.213 

0.214 

0.304 0.210 

0.207 

0.309 

0.310 

1314 

ISO 

.013 

.013 

.037 

.039 

.011 .013 

.009 

.011 

.014 

.039 

D 

0.999 

1.344 

0.979 

1.007 

1.479 1.371 

0.034 

1.187 

0.408 

1.114 

•* FSD 

.891 

.139 

.901 

.334 

.393 .173 

.579 

.330 

.37« 

. 334 

« V D MB 

3.837 

3.079 

2.723 

3.294 

2.374 2.109 

2.339 

1.935 

1141 

1019 

S&r 

9.838 

3.433 

3.401 

3.371 

4.152 3.477 

1357 

1123 

3.799 

3.133 

i«D 

.193 

.071 

.343 

112 

.300 .or 

.133 

.104 

.177 

.099 

« Fan Than 

370 

1030 

310 

930 

480 1000 

390 

070 

430 

930 

« Scattering 

31487 

77904 

31775 

90099 

39133 79311 

33310 

74533 

mu 

48931 

IffldNqi 

0.733 

1.939 

0.549 

0.899 

0.530 1.331 

1.391 

0.853 

1759 

1.804 

terror 

-4.690 

-0.44 

0.305 

4.019 

-31.930 -3.024 

30.839 

9.302 

17.960 

1049 

D 

0.311 

0.213 

0.207 

0.304 

0.209 0.308 

0.309 

0.219 

1331 

1313 

* FSD 

.013 

.015 

.030 

. 028 

.013 .012 

.013 

.019 

.60S 

.035 

P 

0.789 

1.334 

1.948 

1.079 

1.390 1.174 

0.904 

1.353 

1931 

1. 195 

•* FID 

.781 

.174 

.374 

.373 

.344 .190 

.413 

.104 

• SM 

.310 

Dl T* D n* 

2.445 

2.379 

2.077 

1.994 

3.492 1138 

1904 

1004 

1478 

1.910 

Jo. 

3.314 

9.905 

9.739 

9.043 

1943 1399 

3.912 

1399 

3.399 

9.113 

P IB 

.373 

.009 

.213 

.133 

.127 .093 

.305 

.114 

.in 

• 

BBteTtate 

390 

1830 

390 

930 

830 990 

379 

948 

479 

978 

7 gnnSlHteg 

31097 

77907 

31718 

90917 

31953 75913 

11308 

73181 

39989 

88881 

* riBiiinij 

91944 

1.389 

0.999 

0.719 

1.149 1.909 

1199 

0.818 

178V 

2.88V 

%Brrnr 

9.493 

-5.79 

•9.90 

10.71 

-11 734 3.198 

-1011 

133V 

1904 

1888 


® POOR ottatbp 
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However, it should be noted that the dose calculated for the first mode - 

first iteration, which is done with no biasing is a particularly good result. 

Other calculations, which were identical except for the starting random 

number, are not as close to the DOT -HI value as the one shown. The 

results of Table n are based on 200 neutrons for each of the 20 batches 

in each mode calculation. In all cases, except the transport kernel biasing 

-19 

case, the initial total dose is 3. 312 X 10 Rads per source neutron. This 

gives an error of 2. 82 percent. The transport kernel biasing case required 

a different number of random number, so that the initial total dose has an 

error of -13. 47 percent for the same number of neutrons. An even more 

-19 

erroneous result of 3. 976 X 10 Rads per source neutron was computed 
for a run in which 60 batches were used. This represents an error of 
-16.67 percent, even though three times as many neutron histories were 
executed. Thus, the apparent increase in error from the initial to the 
final estimate of the total dose is not necessarily real. In feet, there does 
appear to be a slight decrease in error (and increase in efficiency), but 
further calculations and analyses arc required to substantiate this conclusion. 

One problem with the iterative forward-adjoint method is an increase 
in the error (and decrease in efficiency) in the second mode of the first 
iteration and the first mode of the second iteration. The apparent reason 
for this effect is that the importance function is poorly defined due to the 
small number of histories contributing to the estimate of the importance 
function. Analysis of the importance functions generated by IFAM indicate 
substantial differences in the values of the importance function in phase 
cells which should have identical estimates (e. g. the source region 
importance function for a given energy group in symmetric angular bins 
such »s >, 3, 4, and 5). Figure 10 (A) and (B) seems to indicate that 
this problem is less severe in the latter iterations which have data from 
a large number of histories with which to generate the importance function. 
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Another interesting feature of the iterative forward-adjoint method 
is that the dose estimates tend to oscillate about the correct value. This 
behavior is obvious in Figure 10. Unfortunately, the errors are not 
consistently smaller, but this behavior can be explained on the basis of 
the statistical deviations from the correct dose value. Further analysis 
is required to assure that the long term trend is to smaller error values. 

It is possible that the importance function may cause fluctuations that 
hinder the convergence to a very accurate results rather than assist that 
convergence, as predicted by the "perfect game" in Chapter 2. This effect 
may be caused by "bad actor” histories, which produce a much larger 
importance in a given phase cell than that particular phase cell actually 
has. These "bad actors" have been noted in detailed studies of the energy, 
angle, and region dependent importance functions generated by IFAM. 

Some of these importance function values can be an order of magnitude 
above the expected value. 

Figure 10 (C) and (D) are plots of the total dose estimate in which 
the contribution from collisions in the exclusion volume has been removed. 
While the resulting dose still oscillates, the amount of oscillation is much 
less than for the total dose. This is especially true for the case where only 
the transport kernel biasing was used and where all of the methods described 
in Chapter 4 were in effect. As part (D) of Figure 10 illustrates, the wild 
fluctuations in total dose were largely due to collisions (or lack of collisions) 
In the exclusion volume. Examination of the exclusion volume estimator 
revealed no apparent deficiencies, but further consideration is required. 

The validity of iterative forward-adjoint method as a variance reduction 
technique is at this time questionable. However, due to the fact that toe 
Increase *n ran time aver the time required for MORSE calculations is only 
a few percent, IFAM can be used to great advantage in many problems. 
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One major advantage is that IFAM can be used to compare run time, 
efficiency, or other measures of ''goodness" between the forward and 
adjoint modes of a given problem. Table n illustrates how the adjoint 
mode requires less run time for a given number of histories than the 
forward mode. The percent error for the no biasing case is also 
comparable for the forward and adjoint modes. Therefore, the adjoint 
mode is the better one to use to calculate the air cylinder dose. IFAM 
can be used to aid the researcher in determining the proper mode for a 
given problem. Since very little time is required to set up most problems 
in the adjoint mode, once the forward mode data has been designed, it 
seems advisable to make a preliminary IFAM calculation with both modes 
before selecting the single mode which will be used in the final calculation. 

Another advantage of IFAM is that the importance function generated 
usually agrees quite well with calculations made by more exact methods 
(e. g. ANSIN) once the statistical fluctuations have been removed. Thus, 
the forward and adjoint importance functions from an initial one iteration 
IFAM calculation can be used to suggest subsequent biasing parameters 
for MORSE or IFAM calculations. 
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APPENDIX A 

IFAM Structure and Operating Instructions 

IFAM is a very large and complex piece of computer software. 
Depending on the particular problem, over 100 subprograms may be 
required for a given run. The complexity of IFAM is illustrated by 
Figure A-l, which depicts the interrelationship of the main routine, 
called IFAM, and the various subroutines and functions. Only system 
mathematic library and intrinsic functions have been omitted. Also, 
the subroutines may have more than one version. All cf these sub- 
programs have been included on a FUR PUR tape for the IFAM user. 

If the user wishes to use a specific version, he can easily use the @MAP 
to input a set of source language control statements to specify the 
version. 

A listing of the control statements used to run an IFAM problem 
is given in Figure A -2. Note that the main segment of our overlay 
contains both the main routine, called AAFAM, and the MORSE 
subroutine. There are four first level segments: (1) DA TAIN, 

(2) CRANK, (3) FAMSEG, and (4) NOT USED. DA TAIN is a segment 
that performs all of the input functions except for those few handled by 
MORSE. It is further segmented into the specific random walk problem 
data segment, RWALK, the cross section data segment, XSECIN, and 
the analysis data segment, ANALYSIS. These three segments are 
represented in Figure A-l by subroutines INPUT 1, INPUTS, and INPUTS, 
respectively. This can also be seen in the control statement list by the 
fact that the IN control statement after each SEG statement has these 
element names. IFAM was broken into three input segments so that 
the CRANK segment would dominate the core storage requirements. 

The CRANK segment handles the random walk calculations and the 
o utput . The output can be handled in a separate segment, but the core 
storage savings are small. Since the output routines are called at the 
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end of each mode, it was decided that the core storage savings were 
less important than the extra time for segment interchange. Most of 
the changes made in the Monte Carlo procedure are in the CRANK 
segments. 

The manipulation of the importance information between mode 
calculations is handled by the third segment, FAMSEG. The driver 
subroutine is FAMS, as shown by the element name. (With few exceptions, 
the element name for a given subroutine name are the same. ) The 
fourth segment is made up of subroutines that are referenced in one 
of the above subroutines but are never caUed due to the IFAM problem 
options. Those subroutines include the albedo, fission, and secondary 
gamma rays options as indicated by the IN statement in Figure A -2. 

Thus, this segment is called NOTUSED. By placing subroutines that are 
referenced but never called, core storage is saved for essentially 
loss of computer run time. 

Ao mentioned above, IFAM has some subroutines with multiple 
versions. Most of the original and unchanged subroutines have no 
version identifier, so they can be considered as having die version 
name of "blank. " Those subroutines modified for the IFAM test bed 
were given a version name of IFA . New subroutines written for IFAll 
may have either the "blank" or IFA version name. The CLAS8 control 
statement has been used as shown in Figure A -2 to tell die collector 
that the IFA version of subroutines with multiple version names should 
be included in the collection. A listing of the subroutines stored on the 
FURFUR tape will also show that some subroutines are available with 
DLP version. Subroutines with the DLP version are used in conjunction 
with the IFA versions to reduce the large amount of output, especially 
of the input data. The DLP version is especially useful daring coding 
checkout where the same input data is used several times and a loll 
o u t put each time is unnecessary. 
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The execution of an IFAM job Is illustrated by Figure A -3. The 
stream of control and update cards is representative of that required 
to execute the sar'.ple problem (see Appendix C). The tape which is 
assigned as FAMT contains 3 files, the first of which is a FURPUR 
file and then next two are data files, written with an @DATA statement. 
However, before the FURPUR file can be copied into the temporary 
program file, TPF0, more mass storage must be assigned to TPFS. 

This is accomplished by first releasing it ©FREE TPFjf. ) and then 
assigning a new TPFSf is sufficient FA ST RAND area. 

Next, the source and relocatable elements on the FURPUR file 
are copied into TPF& (If only an absolute element is to be read in 
and executed, the default TPFSf storage is sufficient). This tape read 
is followed by the assignment of a data file DTFF for the forward mode 
input data, which is read from the second file of the FAMT tape. This 
file was written with an ©DATA statement, and consists of card images. 
During the read, columns 52 through 55 in the fourth card are corrected. 
Next, another data file is assigned and the adjoint mode input data is 
read from the third file of FAMT and the sac « correction made to the 
second card. Since the program and data are now on FASTRAND, the 
tape (and tape drive) can be released with a ©FREE. 

The next 5 cards in the control stream are for processing the source 
programs on TPYft. The main routine is updated to reduce the size of 
blank common and the source and relocatable subprograms of CPUTIM 
is deleted. This allows the system routine, C PUTIM/MS FC, to replace 
Hie one on the FURPUR tape. 

After preparing an entry point table, the ©MAP statement is used 
with its control statement? to determine the final executable program. 
FAMSEG Is the name of a source language control element on TPFff 
(and hence on the FURPUR file of FAMT). This element is to be updated 



' by the addition of a library (SYh£s*MSFC#) and the deletion of control 
statement four. The result of these corrections to FAMSEG was 
shown in Figure A-2. A table of the elements in TPF0 is listed as 
a result of the ©PRT, T statement. 

As explained in Chapter 4, data files 16 and 17 are required for 
every IFAM run, and they are assigned to FASTRAND as shown. Now 
the execution of the absolute element FAMSEG generated by the collector 
is initiated. The input card images stored on DTFF and DTFA are 
added to the run stream by the two ©ADD statements. The result of 
this operation will be the IFAM test bed as illustrated in Appendix C. 



APPENDIX B 

IFAM INPUT INSTRUCTIONS 


The input data required br IFAM is essentially the same as 
that required by the Combinatorial Geometry Version of MORSE (see 
Reference 13 ). However, because IFAM requires both forward 
and ajoint mode data for a single run and because some of the MORSE 
options are not implemented inIFAM(e. g. coupled neutron -gamma ray, 
fissions and albedo calculations), this appendix contains a complete set of 
Input instructions for IFAM. Input data which maybe different in IFAM 
than in MORSE has been caveated by placing as asterisk before the 
input symbol in the input data format tables. The change or restriction 
is given in the description for that input symbol or in supplemental 
notes. The input data has been divided into four sections: 


B.l 

Random Walk and Iteration Data 

(B. 1 of Ref. 

13) 

B. 2 

Combinatorial Geometry Data 

(B. 3 of Ref. 

13) 

B.3 

Cross -Section Data 

(B. 4 of Ref. 

13) 

B. 4 

Analysis Data 

(B. 5 of Ref. 

13) 


Each Bisection below contains an introduction to the data format tables 
and defines the subroutines in which the data is read and the differences 
between the forward mode (or F-mode) and adjoint mode (or A -mode) 
input. 

IFAM is so structured that the F-mode data is required to be 
rend initially. Two new data cards are required before the old B. 1 
MORSE input data is rend. These cards are now included an part of the 
new B. 1 F-mode input data. After reading and performing necessary 
calculations with this data, IFAM rends the F-mode date (Mined in B. S, 

B. 3, and B. 4 in order. This data in processed and then stored on temporary 
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data file unit 17. Then IFAM reads the A- mode data, beginning with Card 
Type A of Table B. 1. To enhance the clarity of the output, all title cards 
should specify the input data mode (i.e. F or A). In order to reduce the amount 
of input data required, certain parts of the geometry (B. 2) and cross section (B. 3) 
data is stored ontemporary data file unit 16 during F-mode input ain. ts retrieved 
during the A -mode input instead of reading this redundant data from cards. The 
data handled in this way is explained in Tables B. 2 and B. 3. The input 
and output logical unit numbers have been set to 5 and 6 respectively. 

These units are defined both in MAIN and in Gl, and must be changed in 
both subroutines if other unit numbers are required. 

As explained in Appendix A, it is possible to execute the Com- 
binatorial Geometry Version of MORSE by supplying the proper MAP 
statements to the UNIVAC-1106 Collector under EXEC-Vm. The input 
data defined in this Appendix contains all the necessary information for 
MORSE data preparation. Cards 1 and 2 in the B. 1 sectionapply only to 
IFAM and should not be input for MORSE runs. 

B. 1 Random Walk and Iteration Data 

This section contains all input data for IFAM or MORSE except 
the geometry, cross-section, and analysis data. The first three input 
cards (labelled 1, 2 and 3) are peculiar to IFAM and are not required by 
MORSE. They are read by the MORSE/IFA subroutine for the F-mode 
input only. MORSE then calls INPUT which calls INPUT 1/1 FA. INPUT 1/ 

IFA jandles all of this input for the data defined in this section 

H a description of the input data contains "IFAM:", then any 
information following this symbol pertains only to IFAM input. 
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TABLE B. I - RANDOM WALK AND ITERATION DATA 


Card 

Input 


Type Field/ For mat 

Symbol 

Description 


1 1-5/15 *1MOD I FA M: Initial mode to be executed 

0: F-Mode 
1: A -Mode 

6-10A5 *NITP IFAM: Number of Iterations for this run 

(£ 15) 

11-80/1415 *NPPB(I)| IFAM: NITP order pairs giving the 

1-80/1615 *NBTI(I) ( number of particles per batch (NPPB(I) 

and the number of batches (NBIT) for the 
I-th iteration. 


IFAM: Source region number for the 
F-Mode. 

IFAM: Source region number for the 
A-Mode. 

IFAM: Search length constant for the 
transport kernel routine (NXTCOL) in 
mean free paths. 


3 Intermediate print and mode data write options executed at the end 
of each mode. INTPR(l) - INTPR(12) are for selected labelled 
commons and blank common. INTPR(13) and INTPR(14) control mode 
data writes to data file. Options are: 

0 - Off. Do not write 

1 - On. Write information on appropriate output unit. 

-1 - End. Terminate search for any further output determined 

by INTPR(l) through INTPR(12). Doesn’t affect 
INTPR(13) and INTPR(14). 

I- SA5 INTPR(1) IFAM: Labelled common APOLLO 

6-10/15 INTPR(2) IFAM: Labelled common USER 

II- 15/15 INTPR(S) IFAM: Labelled common OOMLOC 

16 -20 A 5 INTPR(4) IFAM: Labelled common LOC8IG 
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1-5A5 

♦NSRF 

6-10-15 

*NSRA 

11-20/E10.4 

PSC 





TABLE B. 1 - RANDOM WALK AND ITERATION DATA (conk’d) 


Card 

Type 


3 


A 


B 


I 


Input 

Field/Format Symbol Description 


21-25/15 

INTPR(5) 

IFAM: Labelled common PDET 

26-30A5 

INTPR(6) 

IFAM: Labelled common BNKNMC 

31-35/15 

INTPR(7) 

IFAM: Labelled common BANK 

36-40/15 

INTPR(8) 

IFAM: Labelled common RANDOM 

41-45/15 

INTPR(P) 

I FAM: Blank common from 1*1 to 
4*NMTG 

46-50A5 

INTPR(10) 

I FAM: Blank common containing the 
Region Importance (RI) 

51-55A5 

INTPR(ll) 

IFAM: Blank common containing source 
region FAI values 

56-60A5 

INTPR(12) 

IFAM: Labelled common FAM 

61-85A5 

INTPR(13) 

IFAM: Forward Mode output on data 
file 20* 

66-70A5 

INTPR( 14) 

IFAM: Adjoint Mode output on data 
file 20* 



*lncludes ail of blank common and the 
following labelled commons: APOLLO, 
USER, GOMLOC, LOGSIG, PDET, 
BNKNMC, BANK, RANDOM and FAM 

1-80/20A4 

TITLEC; 

Title Card. Any character other than 0 
blank or alphanumeric in Column 1 will 
terminate the job. 

1-5A5 

+NSTRT 

Number of particles per batch 

IFAM: 1NSTRT Is overridden by NPPBfl) 

6-10-15 

NMOST 

Maximum number of particles allowed 


for in the bank(s); may equal NBTRT If 
there is no splitting, fission, and secondary 
generation during execution. If bank slat 
is exceeded by more than 50 due to fission 
or secondary gamma ray generation, the 
job is terminated. 
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TABLE B. 1 - RANDOM WALK AND ITERATION DATA (coat'd) 


Card 

Type 


B 


Input 

Field/ For mat Symbol 


Description 


11-15/15 

16-20/15 

21-25/15 

26-30/15 

31-35/15 

36-40/15 

41-45A5 

46-50/15 

51-55/F.O 

56-60/15 

61-65/15 


♦NITS 

♦NQUIT 

NGPQTN 

NGPQTG 

NMGP 

NMTG 

NCOLTP 

♦IADJM 

AXTIM 

MKDLA 

♦MEDALS 


Number of batches. 

IFAM: Overridden by NBTI(I) 

Number of sets of NITS batches to be run 
without calling subroutine INPUT. 

IFAM: Restricted to 1. 

Number of neutron groups being analysed. 

Number of gamma -ray groups being 
analysed. 

Number of primary particle groups for 
which cross sections are stored; should 
be same as NGP (or the same as NGG 
when NGP=0) on Card XB ready by sub- 
routine XSEC. 

Total number of groups for which cross 
sections are stored; should be same as 
NGP + NGG as read on Card XB read by 
subroutine XSEC. 

Set greater than sero if a collision tape 
is desired; the collision tape is written 
by the user routine BANFR. 

Mode Switch, s0 for F-Moa»: >0 for A -Mode 
IFAM: Overridden by IMOD. 

Maximum clock time in minutes allowed 
for the problem to be on the computer. 

Number of cross-section media; should 
agree with NMED on Card XB ready by 
subroutine XSEC. 

Albedo scattering medium is absolute veins 
of ME DA LB: if 

= 0, no albedo information to be road in, 

<0, albedo only problem — no cross 
sections are to be rend, 

>0, coupled al>edo and transport probl e m 
IFAM: Bet * 0 only. 
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TABLE B. 1 - RANDOM WALK AND ITERATION DATA (Coat’d) 


Card 

Type 


C 


Field/ For mat 

Input 

Symbol 

Description 

1-5/15 

ISOUR 

Source energy group if > 0, if ISOUR « 0, 
SORIN is called for input of Cards El and 
E2. 

6-10/15 

NGPFS 

Number of groups for which the source 
spectrum is to be defined. If ISOUR £ 0, 
NGPFS * 2. 

11-15/15 

*ISBIAS 

No source energy biasing if set less than 
or equal to zero; otherwise, the source 
energy is to be biased, and Cards E2 arc 
required. 

I FAM: Source biasing is handled by code. 
Set ISBIAS > 0. 

16 -20 A 5 

NOTUSD 

An unused variable. 

21-S0/E10. 5 

WTSTRT 

Weight assigned to each source particle. 

31 -40/E 10. 5 

EBOTN 

Lower energy limit of lowest neutron 
group (eV) (group NMGP). 

41-50/E10. 5 

EBOTG 

Lower energy limit of lowest gamma-ray 
group (eV) (group NMTG). 

51-60/B10. 5 

TCUT 

Age in sec at which particles are retired; 
if TCUT=0, no time kill Is performed. 

61 -70/E10. 5 

VELTH 

Velocity of group NMGP when NGPQTN > 0; 
i. e. , thermal -neutron velocity (cm/sec). 

I- 10/E10.4 

II - 20/E 10. 4 
21-30/E 10. 4 

X8TRT ) 
YSTRT> 

zstrt) 

Starting coordinates for source particles 
(Maybe overridden by changes in subroutine 
SOURCE). 

31-40/E10.4 

AGSTRT 

Starting age for source particles (see 
above note). 

41-50/E10. 4 
51 -60/E 10. 4 
61-70/E10.4 

♦UTNP) 
•VINP > 
*WINP> 

Source particle direction cosines; if all 
are zero, isotropic direction are chosen. 
IFAM: Selection based on FAI array. 





TABLE B. 1 - RANDOM WALK AND ITERATION DATA (Coat'd) 


Card 

Type 

Field/Format 

Input 

Symbol 

Description 

El 

1-70/7E10.4 

FS(I) 

NGPFS values of FS, where FS equals the 
unnormalized fraction of source particles 
in each group. 

E2 

1-70/7E10.4 

"BFS(I) 

NGPFS values of BFS(I), where FBS(I) is 
the relative Importance of a source in group 
I. Needed only if HOUR ^ 0 and ISBIAS > 0. 
IFAM: BFS values are not input. 

F 

1-70/7E10.4 

ENER(I) 

NMTG values of the energy (in eV) at the 
upper edge of the energy group boundaries. 

G 

(Omit if NCOLTP on Card B SO, otherwise see Ref. IS). 


1-5/15 

♦NHISTR 

Logical tape number for the first collision 
tape. IFAM: NHISTR cannot be the same 
for F- and A -modes. 


6-10/15 

NfflSMX 

The highest logical number that a collision 
tape may be assigned. 


11-46/3611 

NBIND(J) 

An index to indicate the collision parameters 
to be written on tape (J*l, 38). 


47-59/1311 

NCOLLS(J) 

An index to indicate the types of collisions 
to be put on tape. 

H 

1 

1-24/4X^020 

RANDOM 

Starting random number. 

1 

I 

1-5/15 

NSPLT 

Index Indicating that splitting is allowed 

if >0. 


6- 10 A 5 

NOLL 

Index indicating that Russian roulette is 
allowed if >0. 
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TABLE B. 1 - RANDOM WALK AND ITERATION DATA (Cont'd) 


Card 

Input 


Type Field/Format 

Symbol 

Description 


11-15/15 

NPAST 

Index indicating that exponential transform 
is Invoked if > 0 (subroutine CIREC required). 

16-20A5 

NOLEAK 

Index indicating that non-leakage is 
invoked if > 0. 

21-25A5 

IEBIAS 

Index indicating that energy biasing is 
allowed if > 0. IFAM: If IEBIAS < , 
EPR0B wUl be set to (NGPREG)”* and 
IEBIAS to 1 by INPUTlAFA. 

26-30/15 

MXREG 

Number of regions described by geometry 
input (will be set to one if £ 0). If ENDRUN 
is used, a data array relating media number 
to region numbers must be given in a data 
statement in ENDRUN. 

31-35A5 

MAXGP 

Group number of last group for which 
Russian roulette or splitting or exponential 
transform is to be performed. 


J 


(Omit if NPPLT + N1TLL + NPAST = 0; Repeat the data set given 
below until all groups and regions input; Terminate Card J read 
with negative value of NGP1) . 


1-5A5 

NGP1 

6-10A5 

NDG 

11-15A5 

NGP2 

16-20 A 5 

NRG1 

21-25A5 

NDRG 

26-30A5 

NRG2 


31 -40/E 10. 5 

WTHIH1 

41-50/E10. 5 

WTLOW1 

51-60/E 10. 5 

WTAVE1 


From energy group NGP1 to energy group 
NGP2, inclusive, in steps of NDG and from 
region NRG1 to NRG 2, Inclusive, in steps 
to NDRG, the following weight standards 
and path-stretching parameters are assigned. 
It NGP1 - 0, groups 1 to 1IAXGP will be 
used; if NRG1 = 0, regions 1 to MZREG will 
be used (both in steps of one). Usually 
NDG * 1 and NDRG * 1. 

Weight above which splitting will occur. 

Weight below which Russian roulette is played. 

Weight given those particles ourvtviag 
Russian roulette. 





TABLE B. 1 - RANDOM WALK AND INTERATION DATA (Coat'd) 


Card Input 

Type Field/ For mat Symbol Description 


J 61 -70/E 10. 5 PATH Path-length stretching parameters for 

I use in exponential transform (usually 

0 S PATH < 1). 


K (Omit if IEBIAS on Card I S 0) 

1-70/7B10. 4 *EPROB Values of the relative energy importance 
(IG;NREG) of particles leave a collision or source 
in region NREG for IG*1, NMTG. Input 
for each region (NREG«1, MXREG) must 
begin on a new card. IFAM: EPROB is 
assumed uniform if IEBIAS £ 0. 

(IG = 1, NMTG: NREG = 1, MXREG). 


L 


(IFAM: All variables must be sc 0) 


1-5/15 


6-10A5 
11-15 As 


le-aoAs 


♦NSOIIR Set s: 0 for a fixed source problem; otherwise, 
the source is from fissions generated in 
a previous batch. 

♦MIISTP Index for fission problem, if £ 0 no fissions 
are allowed. 

♦NKCALC The number of the first batch to be included 
in the estimate of k; if s 0 no estimate of k 
is made. 

*NORMF The weight standards and fission weights 
are unchanged if £ 0; otherwise, fission 
weights will be multiplied at the and of 
each batch, by the latest estimate of k and 
the weight standards are multiplied by the 
ratio of fission weights produced in previous 
batch. For time -dependant decaying 
systems, NQRMF should be > 0. 


(Omit if MFI8TP S 0 on Card L) 





TABLE B. 1 - RANDOM WALK AND ITERATION DATA (Cont’d) 


Card Input 

Type Field/ For mat Symbol Description 


M 1-70/7E10.4 *FWLO00 Values of the weights to be assigned to 
| fission neutrons (1 = 1, MXREG). 

IFAM: No fission problems allowed. 


V (Omit if MFiSTP < 0 on Card L) 

1-70/7E10. 4 *FSE(IG, Fraction of fission - Induced source 
IMED) particles in group IG (IG - 1, NMGP) 

and medium IMED (begin a new card for 
each medium where IMED = 1, MEDIA). 
IFAM: No fission problems allowed. 


O (Required input only for coupled neutron-gamma ray problem) 


1-70/7E10. 4 *GWLO(IG, Values of the weight to be assigned to 

NREG) the secondary particles being generated. 

NMGP groups are read for each region 
in a forward problem and NMTG-NMGP 
for an adjoint Input for each region 
must start on a new card. IFAM: Ne 
coupled problems allowed. 
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B. 2 Combinatorial Geometry Data 

The combinatorial geometry data is read by the JOMINAFA 
subroutine, except for the region volumes VNOR(I), which are read by 
the GTVLIN subroutine. JOMIN/IFA is called from INPUT! A FA. Dur- 
ing the initial (F-mode) input, all specified data must be in the input 
stream. However, during the next (A -mode) call to JOMINAFA, card 
types CGB (body data) and CGC (input zone data) are omitted. Since 
the data on the CGB and CGC cards must be identical for both the F- 
and A- modes, this data is temporarily stored on data file unit 16 during 
the F- mode input. It is later read by the GENI subroutine and processed 
for both the F- and A- mode calculations. 

Details of the combinatorial geometry package and its utilisa- 
tion are given in Reference 13. However, the information given below is 
sufficiently complete so that anayone familiar with combinatorial geo- 
metry can prepare input data for most problems. Because the combin- 
atorial geometry package was originally written for the SAM codes (Refs. 

22 and 26) and MORSE originally used the 05R geometry package (Ref. 20), 
so confusion in terminology has occurred. The term zone used below is 
the same as the "region" used in the original combinatorial geometry pack- 
age, whereas a region as used in this document corresponds to "region" 
of the 05R geometry package, just as a media corresponds to the 05R 
"media. " However, regions* and media are defined by combining semes 
for IFAM applications, which differs from the construction in the 05R 
package. The term body has the same meaning as in the original combin- 
atorial geometry package, but it has no counterpart in the C5R package. 

The required combinatorial geometry input data is defined in 
Table B. 2. A summary of the input data required for each body type is 
given in Table B. 3. 

•Note: if ENERUN is used to obtain collision density and track length 
per unit volume estimate of fhience, then a data statement in 
ENDRUN most give a relationship betwe en region and media. 

In this case only one medium may be in a region. 
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TABLE B. 2 - COMBINATORIAL GEOMETRY DATA 


Card 

Input 


Type Field/Format 

Symbol 

Description 


CGA 1-5A5 IV0PT Option which defines the method by which 

region volumes are determined; if 

IV0PT = 0, volumes set equal to i. , 

IV0PT = 1, concentric sphere volumes 
are calculated, 

IV0PT = 2, slab volumes (1-dim. ) are 
calculated (not operational), 

IV0PT = 3, volumes are input by card type 
CGF 

6-10/15 IDBG If IDBG > 0, subroutine PR is called to 

print results of combinatorial geometry 
calculations during execution. Use only 
for debugging, 

21-80/10A6 JTY Alphanumeric title for geometry input 

(columns 21-80). 


CGB (One set of cards is required for each body and for the END card. 

Leave columns 1-10 blank for continuation cards. See Table B. 3 for 
more input information. I FAM- Not required for A -mode input ). 

I- 5/2X, A3 *ITYPE Specifies body type or END to terminate 

reading of body data (for example BOX, 
RPP, ARB, etc. ). Leave blank for 
continuation cards, 

6-10/15 *IALP Body number assigned by user (all input 

body numbers must form a sequence set 
beginning at 1). If left blank, numbers 
are assigned sequentially. Either assign 
all or none of the numbers. Leave blank 
for continuation cards, 

II- T0/6E10. 3 *FPD(I) Real data required for the given boc'y as 

shown in Table B. 3. 
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TABLE B. 2 - COMBINATORIAL GEOMETRY DATA (Cont’d) 


Card 

Input 


Type Field/Format 

Symbol 

Description 


CGC (One set of cards for each input zone. Input zone numbers are 
I assigned sequentially). 


1-5/2X, A3 IALP IALP must be a non-blank for the first 

card of each set of cards defining an 
input zone. If L\LP is blank, this card 
is treated as a continuation of the 
previous zone card. IALP - END denotes 
the end of zone description. 

6-10A5 NAZ Total number of zones that can be entered 

upon leaving any of the bodies used in de- 
fining this input zone (some zones may be 
counted more than once). Leave blank for 
continuation cards for a given zone. (If 
NAZ £ 0 on the first card of the zone card 
set, then it is set to 5). NAZ is used to 
allocate blank common. 

(Alternate UBIASd) and JTY(I) for all bodies defining this input zone. ) 

i HBIASd) Specific the "0R" operator if requir ed for 

the JTY(I) body, 

JTYfl) Body number with the (+) or (-) sign as 
required for the zone description. 


CGD 1-70/1415 MRXZ(1) MRIZfl) is the region number in which 

the "I* 1 *” input zone is contained (I * 1, to 
the number of input zones). Region 
numbers must be sequentially defined from 1. 


CGB 1-70/1415 MMIZ(I) MMIZffl is the medium number in which 

the "1“*" input zone is contained (1*1, 
to the anmber of input zones). Medium 
numbers mt&t be sequent i ally defined 
from 1. 
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TABLE B.2 - COMBINATORIAL GEOMETRY DATA (Cont'd) 


Card 

Input 


Type Field/ For mat 

Symbol 

Description 


CGF (Omit if IVOPT ♦ 3 In card CGA) 

1-70/7E105 VN0R(I) Volume of the "I th " region a - 1 to 

MXREG, the number of regions). 
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INPUT REQUIRED ON OGB CARDS FOR EACH BODY TYPE 


CO CO 

o o 

H CM H 



coca po co coca coco coco 

*3 *3 'S'S 'S’S *3-3 0*3 

»-« CO T1 CO H N H C4 H« 
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TABLE B. 3 (Cont'd) 
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& 

H 


8 

eS 

e 


I 
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B.S Cross Section Data 

The cross section input data required by IFAM are the same 
as lor the MORSE forward and adjoint modes. In order to reduce the 
amount of input data cards for IFAM, the multigroup cross section tables 
(card type XE) are read only during the initial or F-mode input. At 
that time, these tables are stored on data file unit 18 and are automati- 
cally read by READ8G from unit 16 during the A -mode input instead of 
being read from the input data file. In order to impliment this labor saving 
change, the READSG subroutine was revised, producing the READ6GAFA 
version for IFAM. 

To facilitate interpretation of the output cross section data, 
the title card (XA) should indicate whether F- or A -mode output is being 
performed. The data on card XC for the F- and A-mode input need not be 
identical, except for the IDTF switch, since the A-mode cross section 
setup uses the same data as read during the F-mode input and stored on 
data file unit 16. All other card types (XB, XD, XF, XG) require identi- 
cal data for F- and A-mode input. 
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XA 1-80/20AA TITLE Title card for cross sections. This title 

is also written on tape if a processed tape 
generated; therefore, it is suggested that 
the title be definitive. 


XB 1-5/15 NGS The number of primary groups for which 

there are cross sections to be stored. 

Should be same as NMGP input in MORSE. 

6-10/15 NDS Number of primary downs catters for NGP 

(usually NGP). 

11-15A5 NGG Number of secondary groups for which 

there are cross sections to be stored. 

16-20A5 NDSG Number of secondary downs catters for 

NGG (usually NGG). 

21-25A5 I NGP Total number of groups for which cross 

sections are to be input. 

26-30A5 ITBL Table length, i.e. , the number of cross 

sections for each group (usually equal to 
number of downs catters + number of 
upscatters S). 

31 -35 A 5 ISGG Location of within-group scattering cross 

sections (usually equal to number of 
upscatters + 4). 

3 6-40 A 5 NMED Number of media for which cross sections 

are to be stored— should be same as MEIXA 
input in M0R8E. 

41 -45 A 5 NKLEM Number of elements for which cross 

sections are to be read. 

46-50A5 NMDf Number of mixing operations (elements 

times density operations) to be performed 
(must be £l). 

51 -55 A5 NCOEF Number of cosifidsots lor each element, 

including Pq. 
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TABLE B. 4 CROSS SECTION INPUT DATA (Cont’d) 


Card 

Type 

Field/Format 

Input 

Symbol 

Description 

XB 

56-60A5 

NSCT 

Number of discrete angles (usually 




NC ® EF/ Wal K 


61 -65 A 5 

ISTAT 

Flag to store Legendre coefficients if 


greater than zero. Must be > 0 if 
RKLCOL ie used. 


XC 1-5A5 *ERE6G Switch to print the cross sections as 

they are read if 0, if 0 card sequence 
is not checked IFAM: Cross sections 
are also printed if IRD0G * -99. 

6-10/15 ISTR* Switch to print cross sections as they 

are stored if >0. 

11-15/15 IFMU* Switch to print intermediate results of 

M's calculation if >0. 

18-20/15 IM0M* Switch to print moments of angular 

distribution if > 0. 

-85/15 IPRIN* Switch to print angles and probabilities 

it >0. 

26-30/15 I PUN* Switch to print results of bad Legendre 

coefficients if >0. 

3 1-35 A 5 *IDTF + Switch to signal that input format is DTF- 

IV format if > 0; otherwise, ANISN format 
is assumed IFAM: Must be came for F- 
aad A -modes. 

3 6 -40 A 5 IXTAPE Logical tape unit if binary cro s s s e ction 

tape, net equal to 0 if cross sections are 
from cards. If negative, then tbs pro- 
cessed cross sections and other n ecessa ry 
data from a pr evious rua will ba read; in 
this case (IXTAPE <0) no cross sections 
from cards and no mixing cards may be 
Input. The absolute value of IXTAPC is 
the logical taps unit. 
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TABLE B. 4 CROSS SECTION INPUT DATA (Cont’d) 


Card 

Type 

Field/Format 

Input 

Symbol 

Description 

XC 

41-45/15 

JXTAPE 

Logical tape unit of a processed cross- 


section tape to be written. This processed 
cape will contain the title card, die vari- 
ables from common I^CSIG and the per- 
tinent cross sections from blank common. 

46-50A5 I06RT Logical tape unit of a point cross-section 

tape in 06R format. 

51-55/15 IGQPT Last group (MORSE multigroup structure) 

for which the 06R point cross sections 
are to be used (NMGP). 

^Switches are ignored if IXTAPE <0. 


XD (Omit if IXTAPE 0 on card XC) 

1-70/1415 NSIG(I) Element identifiers tor cross-section 

tape. If element identifiers are in the 
same order as element on tape, the 
efficiency of the code is increased due 
to fewer tape rewinds. 


XE (Omit if IXTAPE 0 on card XC) 

Cross sections in A MSN format if 
I DTP 0, otherwise, DTF-IV format. 

Cross sections for INGP groups with 
a table length of ITBL tor NELEM elements 
each with NCOEFF coefficients. 1 7AM: 

Not required for A -mode input. 


XF (NMIX cards are required. Omit if IXTAPE <0 on card XC) 

1-5/15 KM Medium number (media numbers from 

1 to ME HA, see Card B, must appear 
on some XF card). 
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TABLE B.4 CROSS SECTION INPUT DATA (Cont’d) 


Card 


Input 


Type 

Field/Format 

Symbol 

Description 

XP 

6-10/15 

KE 

Element number occurring in medium 


KM (negative value indicates last mixing • 
operation for that medium and at least 
one negative value is required for each 
medium). 

11-20/E10. 5 RH0 Density of element KE in medium KM 

in units of atoms/(barn cm). 


XG (Omit if IOSRTsO) 

1-5/15 NX PM Number of point cross-section sets 

per medium found on an 08R tape, 

= 1, total cross section only, 

* 2, total and scattering cross section, 

- 3, total scattering, and v* fission 
cross section. 
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B. 4 


Analysis Data 

The SAMBO analysis package which was written for use with 
MORSE, into the IFAM test bed. The capabilities of this package were 
sufficient that only minor modifications were required in the analysis sub- 
routines. No changes were required in the input data format or processing 
by the SAMBO package. The input instructions given below for the analysis 
data is read by the SCORIN subroutine for both the forward and adjoint 
modes. A complete set of applicable cards is needed during analysis 
input for both modes. 

Input data for the user written routines, INSG0R, SOURCE 
and ENDRUN, will be inserted into the input data deck (or run stream) 
following the AM cards of the analysis input. 
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TABLE B. 5 ANALYSIS INPUT DATA 


Card Input 

Type Field/Format Symbol Description 



1-80/20AA 

IHOIA) 

Title information describing the analysis 
input data (1=1 , 20). 

1-5/15 

ND 

Number of detectors (set=l if 0). 

6-10A5 

NNE 

Number of primary particle energy bins 
to be used, ( must be NE). 

11-15A5 

NE 

Total number of energy bins (set=0 if 1). 

16-20 A 5 

NT 

Number of time bins for each detector 
(may be negative, in which case/NT 
values are to be read and used for every 
detector) (set=0 if NT 1). 

21-25A5 

NA 

Number of angle bins (set*0 if 1). 

26-30A5 

NRESP 

Number of energy-dependent response 
functions to be used (set=l if 0). 

31-35A5 

NEX 

Number of extra arrays of sise NlfTG 
to be set aside (useful, for example, as 
a place to store an array of group trans- 
fer probabilities for estimator routines. 
If the subroutine ENDRliN which outputs 
flue nee estimates from collision and 
track lengths is used, then this number 
must be at least MXHEG 4 2 (see Card 1 
for MXHEG). 

36-40A5 

NEXVD 

Number of extra arrays to sise ND to 
be set aside (useful, for example, as a 


place to store detector -dependent coun te rs). 
If RELCOL is used, then NKXND oust be 
at least 1 for an event counter. 


(ND cards will be reed, i. e. , 1=1, ND) 

I- 10/El 0.4 XD(1)\ Point detector position. The Astaace 

/ between this point and X8TRT, Y8TBT, 

II- 20/B10. 4 YDCl)} ZSTRT of Card D will be need akn« with 

l the initial ace, AG8TRT, to define the 
2 1-30/E 10. 4 ZD(I)I lower limit of the first tine bia. If other 
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TABLE B. 5 ANALYSIS INPUT DATA 


Card 

Input 


Type 

Field/Format Symbol 

Description 


than point detectors are desired, the point 
positions must still be input and can be 
combined with additional data built in to 
user routines to fully define each detector. 


AD 1-80/20A4 LNK* Title or units for total response for 

all detectors. Will be used in columns 
54 through 133 of the title for the print 
of these arrays. 

(NRESP sets of AE and AF card types will be input, with each AE card 
i followed by the NMTG response values on AF card types for each response). 


AE 1-80/2QA4 LABEL Title or units for each total response 
■ for all detectors. 


AT 1-70/7E10.4 RESP Response function values NMTG values 

will be read in for each total response 
function in order of decreasing energy. 
Each set of NMTG response function 
values will be headed by its AE card 
(I-t, NMTG). 

I 

AG (Omit if NE « 1) 

1-80/20A4 LNK* Units for energy - dependent fluence 

for all detectors. 

AH (Omit if NE s 1) 

1-70/1415 IB 


Energy group numbers defining lower 
limit of analysis energy bins (in order 
of increasing group number). NNE 
in card AB must equal NGPQTN and 
NE must be set to NMTG + NGPQTG 
for a combined problem, or else to 
NGPQTG or NGPQTN. 
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TABLE B. 5 ANALYSIS INPUT DATA (Coat’d) 


Card 

Input 


Type 

Fie id/ For mat Symbol 

Description 


AI (Omit If | NT | si) 

I 1 -80/20X4 LNK* Units for time dependent total response 

for all detectors. 


AJ (Omit if | NT |sl or NE S 1) 

1-80/20A4 LNK* Units for time and energy dependent 

fluence for all detectors. 


AK (Omit if | NT | SI) 

1-7Q/7E10. 4 T(IT, TD) NT values of the upper limits of time 

bine tor each detector in order of in- 
creasing time and detector number. 

The vahiss for each detector mast start 
on a new card. II NT on card AB is 
negative, then only NT values are 
read and used tor every detector. 


AL (Omit if NASI) 

1-80/20A4 LNK* Units for angle and energy dependent 

fluence for all detectors. 

AM (Omit if NA si) 

1-70/tElO. 4 COS (I) Angle bins are defined by the upper 

limits of th e costas of the angle (i.e. , 
the cosine of the lover angle limit). 
Thus, the NAth value mast equal 1. 

♦The Input symbol LNK in n generic variable same for data which In stared 
in blank c omm on and addressed by its location instead of the usual mnemonic 
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TABLE B. 5 ANALYSIS INPUT BATA (Coat'd) 


The inpat cards for user written subroutines are Inserted into the input desk 
after this analysis data. IFAM (and MORSE) calls INSCOR from SCORIN, 
then returns to INPUT3 (may be INPUT or INPUT2 in some MORSE versions), 
which calls USERN. USERN can be used for the input of additional data for 
other user written subroutines. 
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APPENDIX C 
IFAM Sample Problem 

The IFAM sample problem was not chosen because it was the best 
example of the types of problems for which IFAM should be used, but 
because of the availability of experimental and analytical data with 
which to compare results. In addition, important features such as a 
point source and point detector, the need to do biasing to get adequate 
results in a reasonable time, and simplicity of the geometry were taken 
into consideration. These features are present in the point source/ 
point detector configuration shown in Figure C-l. Another attractive 
feature that has been used in testing IFAM is that data for many different 
source/detector separations is available, so the difficulty of the trans- 
port problem can be readily changed. At the same time, the amount 
of effort required to make these changes is minimum. 

The source for the sample point is a point isotropic fission source 
located at (0,0,0). The detector is also isotropic, has a Henderson tissue 
dose response, and is located at (0, 0, Z), where Z is the source/ 
detector separation. Test cases have been run for Z = 150, 800 , 600, 
and 1200 meters. Note that for the adjoint mode, the source is located 
at (0, 0, Z), and the adjoint detector (or forward source) at (0, 0, 0). 

The medium in the volume is uniform air at a density of 1. 11 grams 
per liter. The fission source spectrum, cross sections, and response 
function used are the same as those used for the benchmark calculations 
described in references 28 through 30, which were used to compare 
results. Of course, allowance must be made for the fact that the geometry 
used in this sample problem is finite, while the results of the above 
references are for infinite geometries. 

In order to obtain spatial dependent importance data, 42 importance 
regions were defined as illustrated in Figure C-l. Regions 1 and 2 
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are spheres centered about the source and detector, respectively. In 
the radial direction, the overall cylindrical geometry was divided into 
an inside cylinder and an outside annular volume. In the axial direction, 
five sectors were defined, one below the source, one above the detector, 
and three equally spaced between the source and detector. To consider 
the angular dependence, both the inside cylindrical and outside annular 
region were divided into quadrants for each axial sector. In this 
manner, forty -two importance regions were defined. Energy and 
angular fluxes and importance data wan stored for each region, as 
described in Chapter 4. The input card images, as listed by the 
control stream of Figure A-3, is shown on pages C-4 through C-8. 

Selected pages of the IFAM random walk and analysis output are 
shown on pages C-9 through C- 53. Output from both the forward and 
adjoint modes are given. 
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Figure C-l. Sample Problem Geometry 
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APPENDIX D 


Derivation of the Integral Emergent Adjuncton 
Density Equation 


In Section 2.4 the integral emergent adjuncton density equation (2-117) 
was developed by taking the adjoint of the value equation (2-112), then 
requiring th it adjunctions travel in the direction of the velocity vector 
instead of opposite to it. The MORSE code document (Ref. 11) also derives 
this equation, but the notation which uses integral operators, difficult to 
follow and sufficient editorial errors exists so that the validity of the 
results are suspect on casual inspection. This appendix contains a parallel 
derivation to the one in the MORSE document, but the integral notation has 
been dropped and a different approach taken to obtain the adjuncton density 


equations. The derivation begins with the definition of the two new 

quantities, H (x, to) and G (x, to), which are denoted H (r, w, t) and 

G (r,oo,t), respectively, in the MORSE document. Defining: 

E 


ina 


H g (x,L) 

H g (x,i) 


= 2j(x)X (x, w) 



f R 

J X^(x+R'w)dR* G (x + Rw,w)dR 
0 K 


(D-D 

(D-2) 


* _ 

As pointed out in Reference 11, since X g (x, to) is a flux-like quantity, 
multiplication by the cross section makes (x , to ) a type of event 
density. Note that Equation D-2 differs from the definition given in the 
MORSE document, although the exponential integral term is redefined in 
the MORSE derivation to be consistent with Equation D-2. 


Substitution of the following two equations from Section 2 into 
Equation D-2 and rearranging certain terms produces Equation D-5: 


Rg(x,i) 



(x+ R w)dR’ (£ + R^,w)dR , 

g 


(D-3) 


D-l 




.V f xff*^(x+Ro>; j|w) x __ 

J w,K/-.«--\ X „(x+Rw,b>) dR, 




g- ^ 2f(i + Ri,i) g' 


(D-4) 


g (*+Rw, w ’) 


+ xf (x + Rw)Z 1 ir g (x + Ro>; ojIqj) 
g X^(x+Rw) 


X*(x+Rw,u)dw J dR. 


( D- ' ' 


Since the Xj (x+ Ri) terms cancel and the last term can be multiplied 
and divided by X^ (x + Ru), Equation D-5 can be written as: 

„ -v r®g/-v - f lf(x+FfL)dR'r <t> __ 

H J*.") “j 2n x > e o |r;(x+r«,«) 

» 0 g 

+ £ j— - R(J ?*» ! ") H ,{x+R w) dZ ] dR • (D-6) 

g' lf(i+R*) g 


Comparison of Equations D-2 and D-6 shows that 3 (x,w) 

D 

can be written as: 


G (x,L) - Rj(x,w) +2 J H^x, w)dw . (D-7) 

g g* xf(*> * 


D-2 



The substitution of Equation D-2 into D-7 yields the Fredholm 
integral equation of the second kind for G (x, a>) : 

D 

G (x,a>) = R^(x,<o) +-£ f Ss f 

g B R* "~if (35 jj 1 

R , 

j if (x+R'w)dR' 

e"J ‘ G (x + Rw,w)dRdw\ (D-8) 

The above equation is identical to Equation 2-116, which was derived 
from the value equation. Equations D-6, D-7, and D-8 are the counter- 
parts of Equations 88, 89, and 90 in Appendix A of Reference 11, ex- 

v _ 

cept that the collision operator, C » (r ), has been 

» i 6 , 

defined by placing the (r) - term outside of the summation over g . 

Also, the transport operator, Tg»(r-»>r' , ft'), must be interpreted in the 
sense of the definition given above Equation 88, and not in the sense 
implied earlier by Equations 65 and 66. With these caveats, Equation 
D-8 will also be called the adjoint emergent particle density. However, 
it should be noted that Equation D-8 is not mathematically adjoint to the 
emergent particle density equation, since the kernels are not adjoint. 

Comparison of Equation D -8 with Equation 2-92 reveals that 
the adjoint emergent particle density equation is nearly identical to the 
emergent particle density equation, so that the logic and coding used 
to simulate the emergent particle density equation could be used to 
simulate Equation D-8 with some modification. The adjoint source, 
R^(x, ui), could be input in the same manner as the forward source, 

S (x, E3). The adjoint collision kernel can be calculated from the 
input multigroup cross sections for forward Monte Carlo calculations. 
However, one troublesome aspect of Equation D-8, as well as all 
adjoint equation (2-110, 2-112, 2-116 and D-6) is the transport of 
the adjoint particles in a direction opposite to the direction vector, w' . 
This is due to the fact that the transport step starts are x + RuJ' and 
terminates at x. 


D-3 



One method of overcoming this inconvenience i3 to define a 


new quantity, G (x, w) , where: 

& 

G (x,oj) = G (x, - <lt) . (D-S) 

* o 

Then the integral ’orm of G (x,o») is just: 

C 


f 

8 8 • sf(i) 


g 


..fsfw ^fW«W^ (i+RV _ -. )dW (D 


10 ) 


If the variable of integration, &>', is set to - &>", and the source term, 
, is defined by: 

o 




(D-ll) 


then, 



(D-12) 


The integral over or any other dummy variable representing the 
angular direction, is taken over non-negative definite Increments. 
Since the probability of scattering from - &> to - 5* depends only on 
the angle between - <u and -5"(a s assumed in Section 2. 1), then 



_o |- 5 ) 


(D-1S) 
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Thus, by changing the dummy variable u' to u , Equation D-12 becomes 


G „(x,<j) = ft 


g 


(fa f C Cl- -*'( — V _ ^ , 


K 

,-J 2f (x-R'w') dR • G (x-Rw', w)dRdo . (D-13) 


The above equation has the exact form as the emergent particle den- 
sity equation, and hence can be simulated by the same logic and coding 
without any modification. Care must be taken in the application of 
the G (x,u) simulations since all directions are reversed or in the op- 
posite sense to the normal adjoint functions. Also, the adjoint source 



tion, which must be taken into account for a non-isotropic response 
function. This can be done very easily in the input. 

As discussed in Reference 11 for Equation 93 in Appendix A, 
Equation D-13, the counterpart to Equation 93, can be thought of as 
the transport of psuedo -particles or adjunctons, which travel in the 
"proper" direction. Comparison of these two equations is restricted 
since the definition of the transport operator in Equation 93 is uncertain, 
but it is assumed to be the same as in D-13, except for obvious 
notation changes. In this case, the equations are the same, with 
Reference 11 using the somewhat confusing convention of a different 
symbol for the angular direction variable to represent the change of 
direction from the adjoint emergent particle density equation. Since 
Equation 93 is called the integral emergent adjuncton density, the 
same name will be used for D-13. 
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